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The Berezinskii-Kosteiiitz-Thouless transition in 2D dipolar systems has been studied recently by path inte- 
gral Monte Carlo (PIMC) simulations [A. Filinov et al., PRL 105, 070401 (2010)]. Here, we complement this 
analysis and study temperature-coupling strength dependence of the density (particle-hole) and single-particle 
(SP) excitation spectra both in superfluid and normal phases. The dynamic structure factor, S(q,cj), of the 
longitudinal excitations is rigorously reconstructed with full information on damping. The SP spectral function, 
A(q, lj), is worked out from the one-particle Matsubara Green's function. A stochastic optimization method is 
applied for reconstruction from imaginary times. In the superfluid regime sharp energy resonances are observed 
both in the density and SP excitations. The involved hybridization of both spectra is discussed. In contrast, 
in the normal phase, when there is no coupling, the density modes, beyond acoustic phonons, are significantly 
damped. Our results generalize previous zero temperature analyses based on variational many-body wavefunc- 
tions [F. Mazzanti et al, PRL 102, 110405 (2009), D. Hufnagl et al, PRL 107, 065303 (2011)], where the 
underlying physics of the excitation spectrum and the role of the condensate has not been addressed. 
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I. INTRODUCTION 

Dipolar bosonic systems are of increasing interest for re- 
cent experiments studying the onset of superfluidity in non- 
ideal Bose systems and its connection with correlation and 
quantum degeneracy effects. Examples include dipolar gases, 
as in recent studies of ultra-cold dipolar gases, 1 6 indirect ex- 
citons in coupled quantum wells^ES j n external electric fields 
or exciton-polaritons in quantum wells embedded in optical 
micro-cavitiespJ^ Due to significant experimental achieve- 
ments, a superfluid regime of polar molecules is expected to 
be observed in the near futurePED There is also an increas- 
ing activity on alkali atom gases where the long-long range 
dipolar interactions can be generated by excitations to high 
energy atomic Rydberg states (one electron is excited to a 
very high principal quantum number). The applied moder- 
ate elect ric fie ld can result in a large polarizability and dipole 
momentP^2l All these achievements have a direct impact on 
understanding the properties of Bose-condensates dominated 
by long-range correlation effects. 

From the theoretical side, the static properties of dipolar 
bosons have been quite extensively analyzed in recent years. 
Strongly correlated phases of dipolar bosons confined in a har- 
monic potential have been studied from first principle sim- 
ulations by Lozovik et al.\^ Nho et al.¥^ Pupillo et al.p^ 
Golomedov et al. 25 and Jain et alf^ 2D homogeneous sys- 
tems have been analyzed by Astrakharchik et alW^znd Biich- 
ler et al. 2& with the prediction that above a critical density a 
superfluid dipolar gas undergoes a crystallization transition. 
Properties of a dipole solid have been addressed by Kurbakov 
et al™ The presence of a finite fraction of vacancies and in- 
terstitials, in incommensurate crystal, leads to a superfluid 
response and quasi-equilibrium supersolid phase. In PIMC 
simulations by Filinov et al. 30 the density-dependence of the 
transition temperature from a superfluid to a normal gas phase 
has been analyzed. The composite dipolar bosons, such as 
indirect excitons (pairs of electrons and holes spatially sep- 



arated in semiconductor bilayers) have been recently studied 
by Boning et alW* with focus on the effective exciton-exciton 
interaction and its consequence for the complete phase dia- 
gram of excitonic system. Relevant experimental parameters 
to observe exciton crystallization, in a semiconductor quan- 
tum well, have been predicted. All these results provides a 
reliable estimate of the critical temperature and degeneracy 
parameters required to observe a combined effect of inter- 
particle correlations and Bose statistics in experiment. 

On the other side, theoretical predictions for dynamic prop- 
erties, in particular excitation spectra of density fluctuations, 
are more challenging if performed on a microscopic level. For 
3D dipolar superfluids at low densities the rotonization of the 
excitation spectrum has been predicted from mean-field anal- 
cigar-shaped trap is a promising candidate for ob- 
servation of a roton minimum, where the dipole-dipole inter- 
action is anisotropic, and partially attractive. In a quasi 2D 
trap (a pancake geometry), a roton has a different physical 
nature and is due to a strong in-plane repulsion of similar ori- 
ented dipoles. This regime assumes a high density and goes 
beyond the mean-field predictions. This calls for a more accu- 
rate treatment. The correlated basis functions theory (CBF), 
originally developed for strongly correlated systems like 4 He, 
has been renewed recently by Mazzanti et alP^ to analyze in 
detail the phonon-roton dispersion of an infinite 2D dipolar 
gas at T = 0. The density range has been identified, where the 
CBF dispersion (which include s thr ee-phonon interactions) 
deviates from the Bijl-Feynmarp2l22l an( j Bogolubov predic- 
tions. The upper bound for the phonon-maxon-roton disper- 
sion, derived from the frequency-sum rules, 39 41 is mainly in 
a good agreement. 30 Important extension of the CBF theory 
to 3D geometry and full anisotropy of the dipole-dipole in- 
teraction has been performed by Hufnagl et alW^ This study 
builds an important bridge between the nature of two kinds 
of rotons caused either by attractive or repulsive tail of the 
dipole interaction. Still the current analyses of the excitation 
spectrum lack the important regime of finite temperatures (in- 
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eluding T ~ T c ) and treat excitation damping in an approx- 
imate way! 43 " 45 ! With the present analysis we aim to fill this 
gap and, in addition, discuss the nature of excitations in Bose- 
condensed systems and the role played by a condensate. 

We start the discussion about the nature of excitations from 
the model developed for a one-component neutral BEC. 

The properties of the density fluctuation spectrum have 
been discussed extensively in relation with superfluid liquid 
4 He. In particular, the temperature (T) dependence of S(q, ui) 
from inelastic -neutron scattering experiments was compared 
vs. Glide- Griffi n (GG) model based on the dielectric function 
formalismF^I^The relation with the dynamic structure factor 
is written as 



S(q,u) = - 
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l-V(q)x(q,w)' 
X(q,u) = A(q,uj)G 1 (q,uj)A(q,uj) + X p (q,u) 
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where \ is the irreducible part of dynamic susceptibility x> 
separated into a singular condensate part (related to single- 
particle excitations) and a regular thermal (multiparticle) com- 
ponent, present also above T c . The coupling to the SP spec- 
trum is described by the vertex function A cx ^/n (T) and the 
single-particle Green function G\. The first term in Q van- 
ishes in the normal phase (T > T c ), with the result \ = Xp, 
with Xp being the density response of weakly interacting 
particle-hole excitations, e.g. treated via a Lindhard function. 

The excitation spectrum described by is expected 

to include for small q a zero sound and acoustic phonons, for 
intermediate q a maxon-roton branch, and for large q a com- 
bination of sharp SP resonances and a broad distribution close 
to the recoil energy e q — q 2 /2m due to the multiple scattering 
of weakly interacting particle-hole excitations. 

To fit the experimental observations for 4 He by Eqs.(|T]>-([3jl, 
i.e. a single phonon-maxon dispersion, one suggests that the 
imaginary parts of x an d G\ share a common denominator 
due to the vertex A(g, w)F^As a result the SP and collective 
excitations can not be distinguished in S(q, u) in the super- 
fluid phase. For T > T c the SP spectrum decouples and only 
the thermal part remains (x = Xp)- The dispersion gets broad- 
ened and remains well defined only in the phonon (low-q) part 
of the spectrum. 

In its simple form the dielectric function model failed to 
predict the double roton feature observed in 4 He. Therefore, it 
was necessary to include the mechanism of the quasiparticle- 
decay processes. This has lead to development of theo- 
ries beyond the roton minimum!^«52l \ n his pioneer work, 
Pitaevski53 considered a semi-emperical ansatz for a single- 
particle Green's function 



2uj q T, 12 (q,to) 



(4) 



which includes the unhybridized Feynman-Cohen phonon- 
maxon-roton branchP^ oj„ 



and the self energy, E = 
n^JxiJ, expressed in terms of the two-(quasi)particle re- 
sponse function xi an d me three-point interaction vertex J. 



The self energy accounts for the typical repulsion by anti- 
crossing of two modes, i.e. the interaction between a single 
and a pair of quasiparticle excitations. As a result the SP dis- 
persion, beyond the roton, is renormalized and bends to the 
energy of two roton minimu m. T his concept has been fur- 
ther successfully developed^ P 6 l 59 lby improving both io ¥ q c and 

In particular, Zawadowski et alp] developed the theory of 
a bound two-roton state. Due to hybridization with this state 
the dispersion splits into two branches. The upper branch con- 
sists of heavily damped density excitations with the energy 
close to the recoil energy e q . The low-energy branch is due 
to the SP excitations and saturates slightly below the energy 
of twice the roton minimum. At low T with increasing q its 
weight is transfered to the upper branch. While Zawadowski's 
model provides a good description for low-temperature 4 He 
data by using several fit parameters, it becomes unsatisfac- 
tory at high temperatures. The expected decrease of the in- 
tensity of the lower (SP) branch, by vanishing of the con- 
densate n (T), is not compensated by contribution of atoms 
above a condensate. This is in a contradiction with the sum 
rule, Si,(q) + Sn(q) — S(q), as the static structure factor 
S(q) for 4 He was found to be nearly temperature independent 
for T ~ T c . 49 Further, in some range of g-wavevectors a fit 
to the T-dependent experimental S(q, cj)-data resulted in the 
repulsive interaction between two roton excitations, 60 oppo- 
site to the negative coupling constant suggested in the original 
model. This point, has been resolved by a more elaborated 
theory^using a T-matrix which allowed the coupling constant 
to oscillate with q. 

The dielectric function models interpret the roton and the 
two-roton state as a specific feature of the SP spectrum cou- 
pled in S(q, lu) in a superfluid phase. The standard procedure 
to follow is to fit the involved key parameters (uj q , J, x%) of 
the Green's function Q to be consistent with the low temper- 
ature S(q, w)-data. The T-dependence is neglected and enters 
only via no(T) and the Bose function in Eq. ([TJ). In a similar 
way, the density response function \ p ([3jl is treated also as T- 
independent, and, therefore, is fitted by high temperature data, 
where the singular condensate part does not contribute. 

The weak point of this model is the discrepancy with ex- 
periments near T c . The intensity of the sharp peak (two-roton 
state) assigned with the SP excitations (poles of G\) becomes 
too low as its spectral weight vanishes with n (T). In con- 
trast, a typical experiment on 4 He demonstrates a continuous 
broadening of this peak and nearly conserved integrated spec- 
tral weight. Also a broad multi-excitation background, related 
with x p , extends significantly to low frequencies, as T is in- 
creased, being in contradiction with the model assumption on 
its temperature independence. These discrepancies have been 
reviewed in RefsP^HSSl 

The experiment :) 61 ! 62 ! provided no indication of a well- 
defined mode (SP excitations) suddenly appearing in S(q, to) 
for T < T c due to the first term in Eq. ( 10 1. In contrast, they 



demonstrated that temperature and density dependence of the 
spectral density in the phonon and roton part of the spectrum 
can be fitted with just a single-mode susceptibility. In partic- 
ular, the roton mode demonstrates a rapid attenuation when 
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approaching T c from below and continues as an overdamped 
diffusive mode of zero frequency above T c , The authors con- 
cluded "our results cannot be explained using the GG model, 
unless of course the two components in the GG model hy- 
bridize into one (having one lifetime and excitation energy) 
at all temperatures and pressures, independent of the value of 
no"P^This suggests, that the phonon-roton spectrum must be 
considered as unified branch as was done in the original pa- 
pers by Landau, Feynman and BogolubovP^El However the 
origin of the excitations and their structure can be different. 
By rigorous considerations, Nepomnyashchy 68 has demon- 
strated that a simultaneous presence in the spectrum of both 
zero sound (ZS) and SP branch is not possible. At T < T c the 
ZS is suppressed and replaced by the SP branch defined by 
poles of the SP Green function. The phonons become identi- 
cal with the elementary excitations. The SP branch, however, 
accumulates some properties of the ZS branch, such as linear 
spectrum and independence on uq. This is in a principle con- 
tradiction with the ansatz ([3]) when one uses the substitution 
A oc yJn,a{T). Finally, we note that the Bogolubov type spec- 
trum can be also derived without involving the gauge symme- 
try breaking, i.e. also for a normal phase with the conserved 
gauge symmetry. 69 This implies that the phonon-maxon-roton 
dispersion, observed in a superfiuid, should be also preserved 
in the normal phase, once the damping is small. 

For a 2D dipolar gas, we found that its properties near T c 
qualitatively follow that of 3D 4 He discussed above. From the 
comparison of S(q, uj) at T ~ 0.7T C , we can confirm that the 
intensity of the two-roton "plateau" is unexpectedly too large 
and does not scale with the zero-temperature condensate frac- 
tion once the density is varied. The condensate demonstrates a 
broad variafiorP2] (nn = 0.70 . . . 0.025) and this implies, that 
forAoc y/n (T) in (|}, the two-roton peak should be strongly 
suppressed at high densities, if we assume a roton being solely 
a SP excitation. In our simulations we do not find any evi- 
dence of this behavior. The intensity of the peak is not much 
reduced when approaching T c either. Formally, in 2D geom- 
etry, as considered here, at T ^ the condensate fraction 
should vanish in the thermodynamic limit and no coupling to 
the SP spectrum is expected. Still the two-roton state is ob- 
served for T < T c and even for T > T c at high densities (but 
strongly damped). The finite size effects are negligible for the 
considered system sizes (N ~ 165, 576) and, therefore, do 
not influence the imaginary time correction functions used for 
the reconstruction of the spectral densities (Sec. [Vj. All ob- 
served characteristic spectrum features should remain also in 
the thermodynamic limit (except the limit q —> 0). This brings 
us to a preliminary conclusion that if the two-roton state is 
a feature of SP spectrum the sufficient condition for the hy- 
bridization in S(q, uj) is the presence of at least off-diagonal 
quasi-long range order and a non-zero superfiuid density. 

It is important to mention that, there exists a second view- 
point on roton and different combinations (roton+roton, ro- 
ton+maxon, maxon+maxon) as the intrinsic density excita- 
tions related with the short-range particle correlations tZSHZSI 

To clarify these important issues one needs to indepen- 
dently access S(q, ui) and the SP spectral density A(q, uj) for 
temperatures below and above T c to indentify how the exci- 



tation branches of both spectra are coupled in the superfiuid 
phase and lead to rotonization. This gives motivation for the 
present analysis. The path integral Monte Carlo simulations in 
combination with the stochastic optimization method are used 
to reconstruct collective and single-particle excitations at fi- 
nite temperatures from first principles (in the linear response 
regime). Both should be coupled or form a unified branch 
in the density fluctuation spectrum S(q, uj) once a system has 
a condensate or the Bose broken symmetry. 73 Whether this 
prediction remains also valid for 2D Bose systems, with the 
off-diagonal quasi-long range order and only local Bose con- 
densate, remains an open issue and calls for further analysis 
on the level of microscopic simulations. 

Such analyses in application to 2D dipolar gas will be pre- 
sented in Sec. |VIII| In particular, we found that at weak and 
intermediate coupling (Sec. |VHI A|VIII B"j ) the single-particle 
and the collective (density fluctuations) modes have similar 
dispersion relations in the normal phase (when they are inde- 
pendent). As a result, in the superfiuid phase due to the cou- 
pling (Eq. [3]) both branches become hybridized. In general, 
the dispersion relation uj q observed in S(q, uj) splits into two 
branches, but the high energy (multiparticle) part gets signif- 
icant broadening in the normal phase. At intermediate q (the 
maxon-roton region) the lower branch is shifted to lower en- 
ergies and its damping is enhanced for T > T c . 

For strongly correlated dipolar gas (D = 12.5) our re- 
sults support the viewpoint on the roton as a collective den- 
sity mode. While, the SP spectrum A(q, uj) shows a well pro- 
nounced roton minimum (denoted as £7?), the corresponding 
energies Er ~ 10 and 2Er ~ 20 are larger than those of 
the roton-minimum {Eji ~ 5) and two-roton state (Pitaevskii 
"plateau" 2Er ~ 10) observed in S(q,uj). According to the 
Zawadowski's model, this suggests a strong binding mecha- 
nism which currently lacks experimental confirmation in other 
strongly correlated Bose systems ( 4 He). 

Independent on the presence of the two-roton state, the fre- 
quency sum rules (see Appendix [B]) predict a branch near 
the recoil energy e q . Indeed, the simulations reproduce this 
branch, however, it is strongly damped, as expected from 
weakly interacting density quasiparticles. This multiparticle 
feature is present both below and above T c , In contrast to 
4 He, in a dipolar system at high density (e.g. D = 12.5) 
the Pitaevskii "plateau" appears to be well separated from this 
multiparticle continuum. Similar free-particle like branch is 
observed in the SP spectrum A(q, uj). 

Based on our results we interpret the two-roton state as a 
combination of two density rotons (the poles of the dynamic 
susceptibility) as we can observe it also for T > T c . The 



detailed discussion will be presented in Sec. VIIIB|VIIIC 



In Sec. [n] we introduce the model of a 2D dipolar gas. In 
Sec. Ill physical realizations are listed and experimental pa- 
rameters are compared with our model. In Sec. [IV] the dy- 
namic structure factor S(q, u) is introduced for a system with 
Bose condensate. In Sec. [V] we analyze the structure of the 
imaginary time correlation functions and their relation with 
the SP spectral density A(q,ui) and S(q,u>). In Sec. 



VI 



re- 



construction of A(q, uj) and S(q, uj) using the stochastic op- 
timization procedure is discussed. In Sec. VII we analyze 
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temperature/density dependence of the momentum distribu- 
tion and the spectral densities. In Sec. |VHI| we give our in- 
terpretation of the different excitation branches observed in 
A(q, w), S(q, ui). In Sec. IX we list our conclusions. 



II. MODEL 



The parameters of PIMC simulations are from RefP^ The 
dipole coupling (or density), D ~ n 172 , was varied in the 
range 0.1 < D < 12.5 and temperature 0.7 < T < 3.3. 
With the known D-dependence of the critical temperature®! 
(T c 0.70 ... 1.4) we scan both the superfluid and the normal 
phase. 



Dipolar BEC has been first realized in atomic systems of 
magnetic dipolesP® New exciting results have been recently 
reported for dipolar molecules 14 18 with the major advantage 
of significantly larger electric dipole moments. The mean field 
description of such systems predicts that the dipolar conden- 
sates in 3D geometry become dynamically unstable when the 
dipole interaction dominates the s-wave scattering contact in- 
teraction!22M22lThis effect is a co nsequen ce of the rotonization 
of the Bogolubov-type spectrurrPEEH an( j leads to an insta- 
bility of the time-dependent Gross-Pitaevskii equation. 34 

To avoid this problem and study the behavior of dipolar sys- 
tems well beyond the mean-field regime with a stable con- 
densate and/or a superfluid phase, we consider a quasi-2D 
geometry. The dipole moments are oriented along the direc- 
tion of the strong trap confinement and produce only repulsive 
interaction p 2 /r 3 '. The dipolar forces are assumed to com- 
pletely dominate the contact interaction, therefore, the latter 
is neglected in the present model. We assume 2D spatial ho- 
mogeneity by considering an infinite pancake trap geometry. 
The evaluated spectrum for the in-plane momentum q will be 
valid for the excitation wavelengths, X q = 2n/q, bounded 
by the perpendicular L z and the in-plane L p system size, i.e 
L z < \ < L p . Otherwise, the excitations acquire 3D charac- 
ter. For a quasi-2D model of dipolar condensate, 33 the perpen- 
dicular extension can be estimated by L z ~ l z = \fhfmuj~ z . 

The dipole moment p = ed and particle density n = 1 /a 2 
are free parameters controlled by external fields and in-plane 
trap frequencies. The system Hamiltonian 

N /j2y2 i p 2 

H = -Y.^ n t + 2^e b \r l -r J \3> (5) 

can be made dimensionless by using the length and energy 
units: a = 1/y/n (the mean interparticle distance) and 
Eq = h 2 /ma 2 . In the canonical ensemble, the thermody- 
namic properties of Q are completely defined by the dipole 
coupling D = p 2 /ef,a 3 Eo, temperature T = keT/Eo and 
particle number N. The Bose broken-symmetry condition, 
i.e. (iS?(r)) , (^ + (r)) ^ 0, and evaluation of imaginary-time 
Green's functions suggest to use the grand-canonical ensem- 
ble. In this case, the worm algorithm 75 has been employed by 
fixing a value of the chemical potential /i and system volume 
V instead of N (see note 7 *"). Indeed, existence of the order 
parameter, &o(r) = (^(r)), is related with the fluctuations 
of a particle number in a condensate state. For a weakly inter- 
acting dipolar gas (D = 0.1), we observed significant particle 
number fluctuations, (AiV 2 ) ~ 0.17 (N), with a condensate 
fraction, no ~ 0.7. This can lead to some differences between 
the canonical ensemble (experimental BEC realized in traps) 
and the grand-canonical results. 



III. PHYSICAL REALIZATIONS 

Next we discuss what D values are accessible in current 
experimental realizations. Consider the effective radius of 
dipole-dipole interaction, ad = mp 2 /ebh 2 , corresponding to 
the distance when the dipolar energy reaches the value of the 
zero-point kinetic energy, i.e. p 2 /e&o^ = h 2 /ma 2 d . The intro- 
duced dipolar parameter equals the ratio of two characteristic 
length scales, D = ad/ a. Given a value of a magnetic/electric 
dipole moment (which enters in ad) we can estimate the den- 
sity required for a specific D. The effective dipole interaction 
radius can be expressed as 

a d [k] = 149 ^ 36 m[u]p 2 [Debye 2 ], (6) 

where m is the mass in the unified atomic mass unit, p is the 
dipole moment in the debye. An external electric field aligns 
dipolar particles leading to a repulsive l/r 3 -interaction. Sev- 
eral examples are considered below. 

• Cold bosonic atoms with a permanent magnetic mo- 
ment 1 3 in tight pancake traps. Magnetic dipoles are 
aligned by a magnetic field. A dipolar gas of 52 Cr 
with p = S^b and n ~ 10 n cm~ 2 (a ~ 316 A) 
has dd ~ 24 A. The dipolar interaction will dominate, 
when the s-wave scattering length (a s ~ lOOas) is sup- 
pressed to a s < ad using a high magnetic field Fesh- 
bach resonance. The resulted coupling, D ~ 7.6 • 10~ 2 , 
is close to the analyzed value D — 0.1. The superfluid 
transition temperature is estimated 30 to be T c ~ 12/iK 
(in a 2D geometry). 

Even higher coupling can be achieved in the Bose gas 
of dysprosiurrP ( 164 Dy) with p = 10/ib- The evalua- 
tion for n - 10 9 . . . 10 n cm~ 2 (a - 3160 ... 316 A) 
results in a d - 210.5 A, D ~ 6.6 • 10~ 2 . . . 0.67 and 
T c ~ 0.04 . . . 4/iK. An ultracold dysprosium gas does 
not suffer from chemical reactions, like certain polar 
molecules, 13 and, therefore, is stable at high densities. 

Very recently, a Bose-Einstein condensate of erbium 
( 168 Er) has been reported, 6 being another promising 
candidate for experimental realization of strongly cor- 
related dipolar quantum gases. Similar to dysprosium, 
erbium has a large magnetic moment p = 7/i# and 
atomic mass which enhance the dipole interaction ra- 
dius ([6]). A large dipole moment also allows to reach a 
high efficiency both in a magneto-optical trapping and 
subsequent evaporative laser cooling to temperatures 
T - 200mK. 
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• Polar moleculeP^ (Rb 2 , Cs 2 , 41 K 87 Rb, 40 K 87 Rb) 
transfered to the rovibrational ground state with a per- 
manent electric dipole moment, p = 0.05 . . . 0.6 De- 
byepS aligned by an electric field. For 41 ( 40 )K 87 Rb gas 
with n ~ 10 9 . . . 10 n cm~ 2 and p = 0.57 Debye, we 
estimate a d ~ 6118 A, D ~ 1.9 . . . 19, T C (D = 1.9) ~ 
0.053/zK and T C {D = 15) - 3.4/zK. High densities 
and D-values can be achieved by localizing Feshbach 
molecules in an optical lattice to suppress their colli- 
sions and then quenching to a ground state via stimu- 
lated Raman adiabatic passage. 19 Besides molecule col- 
lisions, an upper bound for quantum degeneracy and the 
dipole coupling D is limited by chemical reactions and 
associated molecule loss rates. 13 

• Rydberg-dressed atoms excited to high principal quan- 
tum numbers k and confined in two dimensions! 21 ! 22 ! 
A large dipole moment can be reached in the linear 
Stark regime with the result p = (3/2)ea#fc(fc — 1) ~ 
1450 Debye for k — 20. We consider an ensemble of 
ground state atoms excited to high Rydberg states with 
the single-atom excitation probability P = (Sl/A) 2 ~ 
0.01 (where f2 and A are a Rabi-frequency and a de- 
tuning laser frequency). The effect of a local dipole 
blockade,! 7 ^! i.e. a zero probability to excite more than 
one atom within a distance R s , can be included via the 
reduced Hamiltonian of "superatoms"pBE3 with the av- 
erage superatom interaction V(r)\ r> R B = p 2 /N 2 r 3 , 
where N s = nR 2 s is a typical number of atoms (in 2D) 
within a superatom radius R s . The latter depends on 
the properties of the excitation laser (P, fl, A), the gas 
density n and the interaction potential V(r). For esti- 
mation, we treat an excitation probability P(n, R s ) of a 
many-body system (modified by correlation effects and 
defined by P(n, R S )N S ~ 1) as that of an uncorrelated 
gas, i.e. P(n,R s ) ~ P ~ 0.01. For a gas of ru- 
bidium Rydberg atoms ( 87 Rb) with the number density 
n ~ 10 9 cm~ 2 we estimate N s ~ 100 (R s ~ 3.1/xm), 
but only a half of the superatoms are interacting as 
they perform Rabi oscillations between their ground 
and excited states. This results in the superatom spac- 
ing a ~ V2R S , and, correspondingly, a d ~ 272/im 
and D = a d /a ~ 61. For low enough temperature a 
Rydberg gas will undergo a crystallization transiti on, as 
the coupling is above the critical value D c ~ I7j 27 l 28 l 3 ° l 
Quantum effects between Rydberg-dressed atoms will 
be relevant for T < 1 nK, due to low densities. 

• Composite bosons, e.g. spatially indirect excitons in 
coupled QWs with the exciton temperature in the range 
0.5 — 4 K. For typical experimental parameters^ 7 ! ( n ~ 
10 10 cm~ 2 ,TO ~ 0.2m e ,£b ~ 10, the inter-well dis- 
tance L ~ 200 A, and p ~ eL ~ 960 Debye) from 
Eq. (JHJl we estimate a d ~ 1500 A and D ~ 1.5. A max- 
imum dipolar coupling (D < 7.5) can be reached be- 
low the excitonic Mott transition density n ~ 1/L 2 ~ 
2.5 • 10 n cm- 2 . 

These estimations show that the systems of ultracold dys- 
prosium (erbium), polar molecules and indirect excitons are 



good candidates for experimental validations of the present 
analyses for 0.1 < D < 12.5. We assume that in an experi- 
ment a long-lived quantum gas is created and, hence, a ther- 
modynamic description remains valid. Further, there should 
be no noticeable heating, e.g. due to breaking of molecular 
bonds or recombination of excitonic states. By this assump- 
tion, experimental parameters which correspond to different 
D-values are listed in Tab. U 

Table I: Relevant parameters for 2D systems of the dysprosium 
184 Dy (Dy), 41 K 87 Rb polar molecules (M) and indirect excitons 
(Ex). The assumed dipole moments are, correspondingly, p Dy = 
10fl B , p M = 0.6 Debye and p Ex = 960 Debye. The dipole cou- 
pling D is expressed in terms of the in-plane density n [cm -2 ]. 
For T < T c a system is in a superfluid phase. 30 The frequency too 
(hujo = Eq) denotes a characteristic energy scale of collective and 
single-particle excitations. 



D 


10 11 


T p y 
£iK 


Dy 

MHz 


n M 
10 s 


rpM 

C 

nK 


KHz 


10 9 


rj-tEx 
C 

K 


GHz 


0.1 


0.023 


0.087 


0.009 


0.027 


0.13 


0.013 


0.044 0.0025 


0.25 


0.5 


0.56 


0.23 


0.22 


0.67 


3.4 


0.33 


1.1 


0.065 


6.33 


1.75 


6.9 


28 


2.7 


8.2 


43 


4.1 


13.4 


0.82 


77 


7.5 


127 


458 


49 


150 


695 


75 


246 


13 


1424 



IV. SINGLE PARTICLE AND DENSITY RESPONSE 
EXCITATIONS IN A BOSE LIQUID 

Existence of a Bose-Einstein condensate is introduced by 
the Bose broken symmetry condition, i.e. by a finite value of 
the order parameter, (&q (r)) = y/Noe 1 ^, where the phase 
4> can be set to zero for a homogeneous condensate. As intro- 
duced by Beliaev 80 the number density operator in momentum 
space can be decomposed as 

P(q) = o.QCi q + at q a + ^2 a+a k+q . (7) 
fc#o 

The first two terms describe the scattering process which cre- 
ates or destroys an atom with zero momentum. The last term 
includes all thermal atoms outside the condensate (k =^ 0). In 
the presence of a BEC the density operator is factorized 

Kq) = \/N^A q + p qi A q = a q + at q , (8) 

which directly demonstrates the coupling of the single- 
particle operator, A q , and the operator of particle-hole exci- 
tations of thermal atoms, p q . The result of this coupling in the 
density fluctuation spectrum was written down by Hugenholtz 
and Pines.^The dynamic structure factor S(q, oj) describing 
propagation of the particle-hole excitations 

i r°° 

S(q,u)=T{( Pq (t) P _ q (0))},T=— J dte^ (9) 

decouples into three parts 

S(q, u) = S A (q,u) + Si(q,u) + S(q,u), (10) 
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where 

S A (q,iu) = T{N (A q (t)A_ q (0))}, (11) 

5i(q, u) = T{ v^Vo[ (p,(t)i- 9 ) + (i,(t)p_ g )]}, (12) 

S( Q ,w)=^{<p,(t)p_,(0)>}. (13) 

Correspondingly, 5a describes the density fluctuations involv- 
ing condensate atoms, Si is a result of the condensate-induced 
intermixing of single-particle and density excitations, and S 
defines the dynamic response of thermal atoms. 

In the long wavelength limit and T = Gavoret and Noz- 
ieresSO showed explicitly that the acoustic phonons are the 
poles of both Sa(<7,w) and S(q,u>). At finite temperatures 
these analyses have been extended by Griffin and Glided 
based on the dielectric function formalism. As a Bose liq- 
uid is cooled down to a critical temperature and a Bose con- 
densate is formed (Nq 7^ 0), the single-particle excitations, 
i.e. the poles of 5a (<7, w), S et a finite weight in the density 
response spectrum. Independent on the representation ([TT)- 
(13 I, as discussed in Sec. [i] the microscopic calculations^ 8 
predict that 5a {q,ui) and S(q,ui) have the same poles in the 
presence of the condensate and that the weight of the SP ex- 
citations should be independent on the absolute value of the 
condensate fraction. 

Hence, the theory allows to interpret the sharp component 
observed beyond the maxon region in superfluid helium as a 
single-particle excitation branch, which vanishes in the nor- 
mal phase with the condensate Nq. This in turn allows to 
indirectly probe the single-particle spectrum of 4 He atoms 
by S(q,uj), as the latter can be measured experimentally by 
neutron-scattering. The development of this concept was 
stimulated by a set of high-resolution measurements!^ 

Our goal is to verify this interpretation and identify the SP 
excitations in S(q, uj). The reconstructed spectra contain pos- 
sible errors coming from the statistical noise in the evaluated 
imaginary time correlation functions and the stochastic recon- 
struction procedure discussed in Sec.|VI| 



V. QUANTUM CORRELATION FUNCTIONS 

In linear response theory an external field produces a weak 
perturbation of a system in thermodynamic equilibrium. In 
this regime dynamical properties can be evaluated via time- 
correlation functions (TCF) of the corresponding dynamical 
operators. For classical many-body systems, TCF can be effi- 
ciently computed in the canonical ensemble by molecular dy- 
namics (MD) simulations. A brute-force solution of classical 
equations of motion nowadays is possible for large ensembles 
(N ~ 10 6 ) of particles. This allows to make accurate extrap- 
olation to the thermodynamic limit, study in detail topologi- 
cal phase transitions^! (e.g. BKT in 2D) involving large-scale 
density fluctuations and divergence of the correlation length 
near a critical point. 

In contrast, for quantum systems first principle time- 
dependent simulations can be done only for few-particle 
atomic and molecular systems, e.g. by time-dependent 
multi-configurational Hartree Fock and CI methods fSl For 




\n-r 2 \/a 2 



Figure 1: Matsubara Green's function (\I/(r2, t2)^ + (ri, ti)} (in 
units a 2 ) for a set of D-values. Simulation parameters (D, fx, V) are 
specified in Tab.|ll| The temperature (T = 1.0) is below to the BKT 
transition. The superfluid fraction is p s /p > 0.80. Due to spatial 
homogeneity we plot the dependence on the relative spatial distance 
\r-2 — ri and relative imaginary time r = t% — t\. 



condensed matter systems the continuous-time Monte Carlo 
method proved to be an efficient and powerful technique.^ 
While predictions from time-dependent methods for equilib- 
rium can be accurately checked by quantum Monte Carlo 
(QMC) methods for ground state and finite temperatures, the 
resolved dynamic correlations are more difficult to control for 
self-consistency. In particular for Bose-condensed systems, 
the importance to accurately treat dynamic correlations be- 
tween thermal and condensate atoms has been discussed re- 
cently by Griffin, Nikuni and Zaremba. 85 

Here, we follow the approch based on quantum correlation 
functions evaluated via path integral Monte Carlo. 

The correlation function of single-particle operators in 
Eq. ( 10 1 (for real time t > 0) contains the contributions of 



two normal and two anomalous Green's functions 

(A q (t)A- q (o)) = (14) 
(a q (t)a+) + (at,(f)a_„) + (fi,(t)a_,) + (fit,(t)a+) 

= -iGu(q,t) - iG 22 (q,t) - iG 12 (q,t) - iG 2 i(q,t). 

The poles of G nm are given by zeros of the common de- 
nominator as follows from the formal solution o f the Dyson- 
Beliaev equation in terms of the self-energies! 73 * 86 ! Hence, 
characteristic single-particle excitations can be extracted from 
one (normal) Green's function Gn 

A(q,u) = A u (q,u) = F{(a q (t)a+)} , (15) 

as other spectral densities A nm share the common poles. 
Therefore, we can concentrate on comparison between 
A(q,Lu) and S(q,u>), and their possible hybridization in a 
Bose condensed phase. 

The general problem with the time-dependent ensemble av- 
erages using stochastic methods, like QMC, is the evaluation 
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\ri - rs|/o 2 



Figure 2: The same as in Fig.[T]but for the density-density correlation 
function (p(r2, t2)p(ri,ti)} (in units a 4 ). 



of the high-dimensional integrals from rapidly oscillating ex- 
ponentials in real-time propagators. This well known "dynam- 
ical sign problem'^ is the main obstacle for first principle 
evaluation of TCF in Eq. ( fl"5] > in the framework of the particle 
orbital-free continuous phase-space methods. 

Another possibility to compute the spectral function is the 
analytical continuation of an imaginary time correlation func- 
tion to real frequencies. 8 ^ The imaginary-time formulation of 
QMC (e.g. DMC, PIMC, GFPIMC) is the best suited tool. In 
particular, the Matsubara Green's function 



G 1 (r 2 ,r 1 ,T) = i'R 



(16) 



can be efficiently evaluated via PIMC in the grand canoni- 
cal ensemble on a set of imaginary time points r G (0, /? = 
\/ksT]. Below we present several examples. 

Fig.[T](Fig.[2]i shows the changes in the structure of the Mat- 
subara Green's function (density correlation function) with 
the increase of the dipole coupling D on the spatial and imag- 
inary time axis. The simulation parameters correspond to a 
superfluid (the superfluid fraction is p s /p > 0.80). In both 
cases, we observe that an enhancement of the inter-particle 
interaction leads to a more complicated structure with oscilla- 
tions. These oscillations decay both in space and time, which 
is a general feature related to damping of one-particle or den- 
sity (collective) excitations of a specific wave-length and en- 
ergy. As we will see later, the interaction also leads to addi- 
tional high-energy excitation branches which are absent in the 
weakly coupled regime, e.g. D = 0.1. 

The behavior of the spatially isotropic function Gi(r, t), 
with r = \r-2 — r\\, can be easily understood in two lim- 
iting cases. As r — > and r — > 0, the many-body ef- 
fects have only a small influence on the one-particle propa- 
gator, which is then close to the free-particle density matrix, 
e -C'-n) 2 /2>i whh A 2 = tf T j m Indeed in Fig. [y near the 
origin we observe a similar Gaussian-shaped peak irrespec- 
tive on the coupling D. After few oscillations in real space it 



Table II: The D and T-dependence of the average particle number 
(TV) and the zero momentum occupation (N q =o) in the volume Vi = 
165. Simulation parameters: Ti = 1.0, T2 = 0.714 and chemical 
potential p,\ . The critical temperature T c is from Ref. 30 The statistical 
error is given in the brakets or if not specified is in the last significant 
digit. 



D 


Mi 




(JV) Ti {N) T2 


(N t=0 ) Tl (N q = ) T2 


0.1 


4.7 


1.30 


163.85(2) 164.82(2) 


107.4(1) 115.4(2) 


0.5 


12.55 


1.35 


163.175(6) 163.308(4) 


77.24(2) 79.52(4) 


1.75 


32.85 


1.39 


169.379(3) 169.403(2) 


44.35(1) 44.95(1) 


7.5 


105.0 


1.22 


165.063(4) 165.075(4) 


9.70(1) 9.885(6) 


12.5 


163.53 


1.01 


164.469(4) 164.492(4) 


3.605 3.885(5) 



Table III: The same as in Tab. [51 for the system volume V2 = 576. 



D 


Mi 


T c 




<^,=o) Tl 


0.1 


4.8 


1.30 


577.85(7) 


342.2(1) 


0.5 


12.70 


1.35 


567.72(6) 


239.3 


1.75 


31.40 


1.39 


563.21(3) 


132.8 


7.5 


103.20 


1.22 


558.06(2) 


27.06 


12.5 


165.53 


1.01 


566.4(5) 


9.048 



evolves into a flat distribution. Another recognizable behav- 
ior is recovered as r — > /3. The Green's function Gi(r,f3) 
becomes the one-particle density with two well-known limits, 
i.e. G\{r = 0, /3) = n, where n = (N) /V is the average 
particle density, and G±(r = L/2,0) being an upper bound 
for the condensate density tiq = (N q= o) /V which decays 
with a power law with system size (in 2D and T 7^ 0). Both 
values depend on fi, V, T and can be read out from Tab. [nj 
The listed values explain the observed off-set in the region 
where G\ (r, r) is almost flat. For temperatures above the crit- 
ical, T c , the off-diagonal quasi-long range order is lost and 
Gi(r,r) decays exponentially with r. With our choice^ of 
the chemical potential fJ,(D) (see Tab. |n| the particle density 
n and the Green's function G\ (0, 0) = n take a similar value, 
i.e. n ~ 1, independent on D. 

The behavior of the density correlator (p(r2, t2)p(ri, ii)) 
(Fig. |2]> has one simple limit. At r = it is a superposition of 
a (5-function (at r = 0) and of the pair distribution function. 
At finite but small t the <5-function turns into the free-particle 
like density matrix, which at strong coupling (large D) gets 
more localized due to the inter-particle interactions. 

Next we consider the reconstruction problem of spectral 
density from the imaginary-time. 



VI. STOCHASTIC OPTIMIZATION METHOD 

In this section we present the stochastic optimization 
method for reconstruction of the spectral densities. The 
method is free of difficulties involved in the analytic continu- 
ation of imaginary time correlation functions. As application, 
the dispersion relations for a 2D dipolar Bose system will be 
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presented in Sec VII 

We start from the general relation between the spectral den- 
sity and the single-particle/density-density propagator Fourier 
transformed to the momentum-space 



G 1 (q,r) = (a q (r)a+(0)) = 



2tt 1 - e-? u 



(17) 



G 2 (q,r) = (p 9 (r)P- 9 {0)) = J du e~™ S(q,u). 

(18) 

The spectral densities satisfy two normalization conditions 

dui 



2tt 



A(q,u) = 1, 



dwS(q,u) = S(q). (19) 



The inversion of equations similar to ( 17 1-( 18 i is known to be 



an ill-conditioned problem and results in a non-uniqueness of 
solution. By Monte Carlo simulations the reconstructed spec- 
tral densities get affected by the statistical errors present in 
G„(g,r). The standard tool used to partially overcome this 
problem is the Maximum Entropy (ME) method!^! However, 
the reconstructed spectral densities for Bose liquids appear to 
be too smooth and the important information on a sharp 8- 
like quasi-particle feature present in the excitation spectra is 
typically lost. 88 Recently, the Stochastic Optimization (SO) 
method has been introduced by Mishchenko et a/., 89 which 
allows to overcome this difficulty. Its main advantages com- 
pared to the standard regularization methods, like the ME, are 
the continuous parametrization in frequency space, instead of 
a predefined mesh, and non-suppression of high derivatives of 
the spectral function, performed by the regularization meth- 
ods. As a result, sharp peaks and edges are not lost during the 
reconstruction. In its core, with the SO one solves by stochas- 
tic sampling the minimization problem of the least deviation 

D n [G„] =X)I 1 ~ G n {q,n)/G{q,n)\w{q,n), (20) 
n 

where 6G(q, Tj) is the statistical error of G(q, ri) at the imag- 
inary time T{ (At = r i+1 — Tj, with i = 1,...N and 
TV At = ($), the weight factor w(q, Tj) is chosen in the form 

w(q,n) = (N/N T )w(q,Ti), N T = W (<1, n), (21) 



Ti 

w(q,Ti) = min[10, \G(q,n) / '8G{q,n)\], 
w(q,n) = max[l, \G(q, n)/5G(q, n)\], 



(22) 
(23) 



and G n is generated from 

G n {q,r) = I e-™S n {q,u)duj, (24) 



with S n (q,uj) being a trial spectral density parameterized into 
some basis set. The deviation D n [G n ] is optimized by a ran- 
dom sequence of updates which can change both the param- 
eters of the basis set and its size. By increasing n (the num- 
ber of independent solutions) we evaluate the corresponding 



variance ajj = (D%) (with the zero mean). In the end, we 
select only "good" solutions from the whole sample which 
satisfy the condition D n < D m [ n with D^ m = 1.5 ct,d„- 
The final solution is constructed from their linear combina- 
tion (100 < M < 400) 



1 M 

S est {q,u,M) = — ^2S n {q,uj), 



(25) 



71=1 



to take advantage of a self-averaging of the noise. 

For the parametrization of trial G n we use a set of rectan- 
gleP 



{-Pi}i=i,jv = {hi,Wi,Ci}i—i t j^-, 

N 



hi{to) 



7^0, u £ [cj - Wi/2,Ci + Wi/2], 
0, otherwise. 



(26) 
(27) 

(28) 



with height h, width w and center c being the optimization 
parameters for a fixed value of the g-vector. The basis size 
TV was also varied during the optimization in the range 80 < 
iV< 400. 

For the reconstruction of the dynamic structure factor 
S(q, ui) we have used only positive frequencies, < uj < 
coco, with the cut-off frequency uco ~ 400, by taking into 
account explicitly the relation between negative and positive 



energy transfers, i.e. S(—q, —uj) 
suits in 



-0U 



S(q, uj). This re- 



Gn(q,r) 



J S{q,Uj) (l + e-^) duj. 



(29) 



In the basis of the rectangular functions the trial imaginary 
time density correlation function takes the form 



G 2 (q,r) 



2 I> [ E l r Cit 

i=l \t=T,P—T 



sinh 



Wit 



(30) 



which is symmetric with respect to the mid-point t = (3/2 
Therefore, it is sufficient to evaluate G 2 (q, r) (see Eq. ( [T8] 
for the imaginary times r £ [0, (3/2]. The normalization con- 
dition ( fT9] l results in the constraint 



G 2 (q,T = 0) = G 2 (q,/3) = S(q), 



N 



^ hiWi +S' = S(q), 



i=l 



TV 



S' = 2^e-^sinh^f, 



(31) 
(32) 

(33) 



where the factor S' corresponds to the occupation of the neg- 
ative frequencies with the meaning of the energy transfer huj 
from a system to a scattered probe particle, i.e. from the ex- 
citations existing in the system. Typically, this contribution is 
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important in the g-region of the acoustic phonons and rotons, 
and one can start the optimization first by neglecting S' and 
find an optimized solution Si(q, uj). Then the normalization 
of a second solution is corrected by 



JV 



J2 h 



(n) (n) 



= S(q)-S'(S est (q,uj,M)), 



(34) 



i.e. the correction S' is evaluated based on the spectral den- 
sities obtained in the previous iterations (M > 1). With n 
increasing the on fly estimation of both S est (q, w, M) and 5" 
is improved. This results in a fast convergence of S' within 
few (n ~ 5) iterations. 

The reconstruction of the spectral density A(q,ui) of the 
one-particle Green's function ( [17} gets a bit more involved. 
First, there are no simple relations between the densities at 
positive (A > (q,ui)\ u >o ^ 0, A>(q, w)| w <o = 0) and neg- 
ative (A < (q, w)| w < ^ 0,A < (q,ui)\ UJ>0 = 0) frequencies. 
They should be worked out independently in the frequency 
range, — ujqo < w < ^co We found that using ujqo ~ 600 is 
sufficient to fit a fast drop of G\ (q, r) near r = 0, see Fig. [5] 
Second, Gi(q,r) is not symmetric relative to r = (3/2 and 
should be evaluated on the whole interval, < r < f3. Third, 
the Laplace transform in Eq. ( [17} contains the additional Bose 
factor which leads to the rigorous resurtp2U<(g,w) < 0. 

A simple solution which allows to use the same SO proce- 
dure as for S(q, uj) is to introduce two spectral densities which 
are both positive and defined for w > 



A > (q,u;)=A > (q,u;)/(l-e 
A < (q,oj)=A < (q,u } )/(e-^ 



~ Pu ) > 0, (35) 
- 1) > 0. (36) 



This results in the following decomposition 

2tt 



Gi(q,r) 







-TW 



A > ( g , W ) + e-^- T ^ < (g J a;) 

(37) 



Using Eqs. ( 17 1,( 19 1 we end up with two normalization con- 
strains 



^A > (q,uj) = G 1 (q,0) 



2i 



-^i<(<7^), (38) 



pA < (q,u) = G 1 (q,!3)- f ^ '£>(<,, a,), (39) 





G 1 (q,0) = l + G 1 {q,p). 



(40) 



At low temperature and/or high excitation energies, j3ui 1, 



the integral terms in Eqs. (38i,(39i can be neglected. Then 



both spectral densities become independent with the normal- 
ization given by the momentum distribution 



■\{m\4\n)\ 2 = 



duj A(q, uj) 
2tt - 1 



The latter can be directly evaluated via Fourier transform of 
the one-particle density matrix 



n(q) = (a q a+) = 



[ dre- lqr ^(r,t). 
Jv 



Ar&r' e iq{r '- r \i>{r' ,t)$ + (r,t)), 



(42) 
(43) 



In contrast, at high temperatures, the normalizations ( 38 1-( 39 1 



are mutually dependent and can be treated iteratively, simi- 



larly to Eq. (34 1, i.e 



]T h$w%> = n(q) + 1 - A'(A estt< (q, «, M)), (44) 
i=i 

N< 

J2 h< A w S=n(q)-A'(A e ^ > (q,u J ,M)). (45) 
1=1 

The procedure converges within few iterations. 

For the detailed description of the stochastic optimization 
method and types of the update we refer to Refill We found 
of particular importance to implement the annealing allowing 
to escape from local minimum and minimization of the devia- 
tion ( |20| using parabolic interpolation. During the reconstruc- 
tion, the updates involving changes in two rectangles ( f26| (i.e. 
change a weight of two rectangles, add a new rectangle, re- 
move a rectangle) correctly redistribute the total spectral den- 
sity A(q, u) between Ay (q, oj) and A< (q, uj), even if initially 
they are chosen to be equal. The decay of spectral weight 
A< (q, lu) at large q-vectors should follow that of the momen- 
tum distribution and, therefore, practically vanishes beyond 
the roton region (qa ~ 6). 

To speed up the optimization process and convergence 
of the finite-temperature corrections in the normaliza- 



tions ( 34 1,( |44|) ,(45 i the centers of the basis set rectangles {cj}, 
representing S(q, uj), Ay (q, uj), A< (q, uj), have been initially 
normally distributed around the frequency ui x (q), i.e. the up- 
per bound of the lower excitation branch, see Appendix [B] 
The optimization process was started with JV, iV> , N < ps 30 
basis functions with the initial width, w — w m [ n , i.e. equal 
to the frequency resolution hw m { n /E^ — 0.5. The optimiza- 
tion was performed in several iterations each of 120000 steps. 
In each step one of the following update types has been ran- 
domly chosen: (1) shift of a rectangle, (2) change a rectangle- 
width, (3) change the heights of two rectangles, (4) add a new 
rectangle, (5) remove a rectangle, (6) split a rectangle into 
two, (7) glue two rectangles. For each update involving a 
change in one or couple of parameters {hi}, {wi}, {c,} we 
try to converge to a local minimum by parabolic interpola- 
tion. In the first 60000 steps we perform several annealing 
sequences with a duration 1000 — 5000 steps and temporary 
accept the updates increasing the deviation D n (20 1. Once an 



iteration is finished, the actual deviation D n and the optimized 
solution are saved and checked for the acceptance criterion, 
= G\(q,p). D n < D m i n . If not accepted, we proceed to the next itera- 
tion simultaneously increasing the minimal deviation by some 
(41) factor, D m [ n —> jD m \ n , e.g. 7 = (1.1 + v), with the random 
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Figure 3: (Color online) (a) Fourier transform of the one-particle 
Green function Gi (16) : PIMC data (symbols) and the stochastic op- 
timization fit (solid lines). Vertical dashed error-bars demonstrates 
the statistical error (shown only for two curves). Colors (and num- 
bers) correspond to the chosen g-vectors shown on the left panel, (b) 
Relative deviation < |20| > after optimization. Simulation parameters: 
D = 12.5 and T = 1.25 (see Tab. |n). Temperature is close to T c . 
Superfluid fraction p s /p = 0.38(2). 



ficulty to distinguish individual curves demonstrates the con- 
vergence of the optimization process. All solutions are ob- 
tained from independent (randomly chosen) initial spectral 
densities. The relative error depends on the g-vector and its 
value is mainly determined not by the optimization result but 
by the statistical fluctuations near r = 0/2, when the value of 
the correlation function significantly drops. Still the optimiza- 
tion process converges. A well-behaved decay at short times 
(for t or (3 — r) fixes some of excitation energies u) in the 
spectral density and, hence, partially the decay at large times 
when approaching the mid-point r = j3/2. As follows from 



Eq. ( 20 1 the fitted points contribute to the least deviation D n 



not equally but with the weight factor determined by the sta- 
tistical error. Hence, the middle points with large statistical 
fluctuations are less important for the fit. 

The influence of the statistical noise on the accuracy of the 
reconstruction procedure is discussed in detail in Appendix|A] 



VII. STATIC AND DYNAMIC PROPERTIES 
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Figure 4: (Color online) The same as in Fig. [3] for the density cor- 
relation function G2 (18) . Simulation parameters: D — 12.5 and 
T = 1.0. Superfluid fraction p s /p = 0.81(2). G 2 (q,r) is evalu- 
ated only for r G [0, 13/ 2] due to the symmetry with respect to the 
mid-point r = /3/2, see Eq. 1 29 1. 



number v G [0, 1). The initial value Dmin was chosen based 
on the statistical noise in the quantum correlation functions, 

10~ 5 . 
10- 4 . 



Eqs. ( 17 1,( 18 1. Typically we start from 



10" 4 
10- 3 



and can reach the acceptance error D m i n 
within 5 — 6 iterations, depending on the level of noise in the 
simulation data. 

Figs.[3]and|4]demonstrate two examples of the performance 
of the optimization procedure. The symbols show the PIMC 
data Fourier transformed to the momentum space and demon- 
strate the level of statistical noise. For the density correlation 
function we get a typical error, 6G2 ~ 10~ 4 , which is one 
order of magnitude smaller compared to the single-particle 
Green's function, 8G\ ~ 10~ 3 . The solid line is the result 
of the optimization, i.e. the correlation function ( |24| evalu- 
ated from one of the solutions S n (q,oj) or A n (q,u) which 
enters in the final estimation (25 1. The right panel shows the 
relative error for a set of solutions G n (n = 1, 100). The dif- 



A. Momentum distribution 

Before analyzing of the spectral densities, we now discuss 
the momentum distribution which enters in the normaliza- 
tion |44|-(|45|. These data will characterize the normal and 
superfluid phase of 2D dipolar gas. 

In Fig.[5}t,b,c the low (T < T c ) and high (T > T c ) tempera- 
ture momentum distribution is shown for D = 0.1, 1.75, 12.5. 
These are referenced in the following as weak, intermediate 
and strong coupling, correspondingly. The T-dependence of 
the one-particle density matrix ( |42) i is shown in Fig. [5jJ,e,f. 
For low temperatures (T — 0.714, 1.0) in the superfluid phase 
it demonstrates only minimal changes. Our system has a finite 
volume and satisfies the periodic boundary conditions (PBC), 
therefore, the momentum distribution is evaluated at a discrete 
set of wavevectors, q = 2ttti/L (n G Z), shown by symbols 
in Fig. [5] For comparison, the solid line is the result (shown 
for T = 1.0,3.3) obtained by extension of Eq. ( p2] > to the 
limit V — > 00 with the assumption that G\ (r, (3) decays as a 
power-law (exponent) below (above) T c beyond the simula- 
tion box size L — \/V . The fitting parameters for T < T c 
are given in Tab. [IV] This method agrees with the finite-size 
results (shown by symbols). As the system size L gets larger, 
more discrete values of q will get occupied dwelling on this 
interpolation curve. 

We get a full agreement in the normal phase, when G\ (r, f3) 
decays fast to a small value at r = L/2 (e.g. to ~ 10~ 6 for 
T = 3.3 for D = 12.5), and finite-size effects are of minor 
importance for systems with V > 165. 

Some discrepancies appear in the superfluid phase. In par- 
ticular, for T = 1, D = 0.1 (Fig. |5^) and qa > 5 we observe 
some statistical noise in n q related with the Fourier transform 
of G\(r, 0). In the long wave-length limit (qa < 0.5) the dis- 
crepancies are induced by the interpolation. The density ma- 
trix is influenced by the PBC and deviates from the expected 
power-law decay r^ v ^ near r — L/2. Therefore, the fit with 
the critical exponent v(T) was applied in the range 4 < r < 6 
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Figure 5: (Color online) Temperature dependence of the momentum 
distribution n q (a-c) and one-particle density matrix Gi (r, /?) (d-f) 
in the log-scale for dipole coupling D — 0.1, 1.75, 12.5. In the su- 
perfiuid phase, the decay of Gi(r, f3) is characterized by: 1) a fast 
short-range decay to a condensate fraction no(T); 2) off-diagonal 
quasi-long-range order which depends on temperature and the super- 
fluid density. Simulation parameters are V = 165, chemical poten- 
tial ^(D) — 4.7, 32.85, 163.53, the average particle number (N) ~ 
164. The used wavevectors: qa = 2nn(a/L), (n — 0,1,...). 
The solid lines are obtained by interpolation (see text). The system 
is partially superfluid for T < 1.67 (D = 0.1, 1.75) and T < 1.0 
(D = 12.5). 



Table IV: The D and T-dependence (Ti 
the critical exponent by fitting the long-range 
particle DM, Gi (r) /n ~ 
volume Vi = 165. 



^ T >,andGf /2 



1.0, T 2 = 0.714) of 
asymptotic of the one- 
= Gi(L/2)/nforthe 



D 


o(Ti) o(T 2 ) 


v(Ti) i/(T B ) 


Gf /2 (T0 


Gf /2 (T 2 ) 


0.1 


0.738(4) 0.762(6) 


0.088(3) 0.063(5) 


0.628(4) 


0.678(3) 


0.5 


0.526(4) 0.523(4) 


0.085(4) 0.061(5) 


0.449(3) 


0.467(4) 


1.75 


0.294(2) 0.285(1) 


0.103(3) 0.075(2) 


0.245(2) 


0.250(2) 


7.5 


0.0656(5) 0.066(1) 


0.159(4) 0.155(8) 


0.0495(2) 


0.050(1) 


12.5 


0.0243(5) 0.025(1) 


0.258(5) 0.22(3) 


0.0153(1) 


0.0172(6) 



for Vi = 165 and 4 < r < 10 for Va = 576. This partially 
allows to exclude the effects of the short-range correlations 
and of the finite-size errors on the expected asymptotic decay 
of G\. As predicted by the interpolation curve, in a system 
with off-diagonal quasi-long range order, the momentum dis- 
tribution diverges when approaching zero momentum. For a 
finite system and discrete q we can only see the onset of this 
regime. For V — 165, the smallest wavevector (qa — 0.486) 



Table V: The D and T-dependence of the condensate fraction no (T) 
(for system size Vi = 165 and V2 = 576) at Ti = 1.0 and 
T2 = 0.714 evaluated from Tab. [ P| |Tll] compared with the zero- 
temperature condensate no(O)! 27 ' 77 ^ 



D 


^(1) 


(0.714) 




710 (0) 


0.1 


0.655(5) 


0.700(5) 


0.592(4) 


0.72 


0.5 


0.473 


0.487(1) 


0.422(2) 


0.50 


1.75 


0.262 


0.265 


0.236 


0.28 


7.5 


0.0587 


0.0599 


0.0485 


0.062 


12.5 


0.0219 


0.0236 


0.016 


0.025 



is too large to capture the divergence. The low-momentum 
behavior will be discussed in more detail later. 

The information about condensate or occupation of the 
zero-momentum state is given by the value of the Matsub- 

0,/3) 



ara Green's function (41 



distribution in Eq. (42 



Gi(q — 0,/?). The momentum 
is normalized by the average particle 
number (N) in volume V. Therefore, the number of parti- 
cles at zero momentum (N q= o) also depends on V and (N), 
see Tab. [V] For a macroscopic system this value will diverge 
as V — > 00, in agreement with the interpolation curves in 
Fig.[5^,b,c. On the other hand, at a finite temperature the con- 
densate fraction should vanish in the thermodynamic limit, 40 
limy_j. 00 (N q= o) /V — 0, as the zero momentum state will 
be depleted by thermal fluctuations. A number of occupied 
phonon modes with q = 2nn/L < q c gets larger with in- 
creasing L. 

How the thermal depletion of the low-momentum states 
proceeds can be analyzed for T = 1.0,2.0 and 3.3. For 
D = 0.1 broadening of n q is observed in a wide range of 
momenta. The formation of the condensate feature at small q 
is accompanied by suppression of the high-momentum states. 
For D = 1.75 this is noticeable for qa < 1, and for D = 12.5 
for qa < 0.5. Here, formation of the condensate is due to 
suppression of the states with qa > 6. We conclude that the 
increase of the coupling/density narrows the interval of mo- 
menta where a fast divergence can be observed in the super- 
fluid phase. This, certainly, will complicate an experimen- 
tal detection of a superfluid transition in a strongly correlated 
Bose system based only on the specific features of the momen- 
tum distribution. In particular, for D — 12.5, in a broad range 
of momenta (1 < qa < 6) the distribution is practically tem- 
perature independent (0 < T < 3.3). Here the main depletion 
mechanism (also at T = 0) is a strong interparticle interac- 
tion. Only a small fraction of particles occupies the q = 
state. The depletion out of the condensate is enhanced by the 
presence of low-energy excitations, e.g. rotons. More insight 
should be given by the spectral densities of single-particle ex- 
citations and their dependence on the interaction strength. 

The average particle number in the zero-momentum state 
(N q=0 ) and the condensate fraction n Q (T) for T = 0.714, 1.0 
(T < T c ) are given inTab^ |ll|Vl The PIMC results are com- 
pared with the DMC. 27 77 91 The zero-temperature condensate 
fraction stays the upper bound for the PIMC values, no (T) < 
nrj(0). Both results strongly deviate from the predictions 92 
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based on the perturbation expansion in the gas parameter na 2 . 
As shown in RefsJMIIlthe s-wave scattering model for the in- 
teraction fails at the densities na 2 > 10~ 3 (D > 0.01) signif- 
icantly lower then considered here (D > 0.1). 

The T-dependence of the condensate is more pronounced 
for D = 0.1 and 12.5. In the first case, it is due to occupation 
of the higher-momentum states (see Fig. |5^). In the second 
case, Ti and T 2 are close to the critical temperature of an infi- 
nite system (T c — 1.01, Tab.|H} and thermal fluctuations play 
an important role. 

The condensate fraction is also estimated for a larger sys- 
tem (V2 = 576) to demonstrate the finite size effect. For D = 
1.75,7.5 and 12.5 we find G^ 2 (T X ) = 0.220(5), 0.0438(8) 
and 0.0136(5), correspondingly (G\ = G\jn). Compared 
to n (T) in Tab.|v] the relation, n (T) > C[ /2 {T), always 
holds and connects the reduction of the condensate with the 
boundary value of G\. The difference of both decreases with 
v{T). For D = 0.1, 0.5 and 1.75 (Tab-jlV)) the difference is 



within few percents. For D = 7.5 and 12.5 for a better agree- 
ment T should be lowered to decrease the critical exponent 

~ L/2 

Taking G x as the estimation of the condensate fraction at 
T ^ (when v is small), some predictions can be made for 
experimental systems with a number of bosons beyond direct 
numerical treatment (see Appendix |C]>. 

Next, we consider the divergence of the momentum distri- 
bution as q —> 0. Small system size limits the resolution in this 
important region to q > 2tt/L. To overcome this limitation 
we employ the "1 /p 2 sum rule'^JEH 



+00 

lim -Gn(p,w = 0)= / — 
p->0 J Zir 



du> A(jj, oj) mriQ 



PsP 



(46) 



which holds independently on the coupling strength. In the 
collisionless (T -C T c , p s = ri) and low-frequency (hydrody- 
namic) regime we can apply the ansatz by Gavoret and Noz- 
ieresP 

A(q,u) = 2n{Z(q) + l)6(w - cq) - 2nZ(q)S(u + cq), 
where c is the compressional sound speed. Substitution of 



this ansatz in ( 46 1 yields both the spectral weight Z and the 
T-dependence of the momentum distribution 46 



Z(q) 



mnQC 



N(uj) 



1 



2 Ps q ' w eP" 
q = Z{q) [2N(cq) + 1] . 



(47) 
(48) 



The sound speed c(T) can be evaluated from the particle num- 
ber fluctuations 



(N 
c = 



(JV) S 



(N) z k B TK T \/V, (49) 
[ran k,t\ ~' ~ , (50) 
where is the isothermal compressibility. The sound speed 
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Figure 6: (Color online) Long wavelength asymptotics of the mo- 
mentum distribution n q in the superfluid phase, Tip) = 1.0(0.714) 
vs. D (log-log scale). The dashed lines are the prediction by 
Eq. l |48| > for 7\(2)- Simulation parameters: Vi(2) = 165(576), 
T m = 1.0(0.714). 



Table VI: The D- and T-dependence of the compressional sound 
speed c(T) [aEo/h], isothermal compressibility kt(T) x 10 2 
[h 2 /m] and the superfluid fraction p s (T) = p s (T)/n(T). Tem- 
peratures are T = 1,2 and T = 3.3 [p s (3.3) = 0]. Simulation 
parameters are the same as in Tab. 1 1 1 1 1 1 [ | Second line for each D is 
for the system volume V2 ■ 



D 


c(l) c(2) c(3.3) 


«t(1) kt(2) 




0.1 


2.366(5) 2.213(5) 2.59(2) 
2.411(7) 2.26(6) 2.72(5) 


17.98(8) 23.4(1) 
17.15(10) 22.0(1.1) 


0.95(1) 0.032 
0.98(3) 0.004 


0.5 


4.097(8) 3.863(9) 4.06(1) 
4.16(2) 3.88(2) 4.12(2) 


6.024(24) 7.05(4) 
5.85(4) 7.0(1) 


0.98(1) 0.10 
0.98(2) 0.007 


1.75 


6.795(8) 6.66(3) 6.88(4) 
6.69(3) 6.46(2) 6.67(2) 


2.110(5) 2.22(2) 
2.28(2) 2.48(2) 


1.0 0.17 
1.00(3) 0.005 


7.5 


12.32(2) 12.38(6) 12.40(11) 
12.25(5) 


0.659(2) 0.655(6) 
0.688(5) 


0.95 0.007 
0.95(2) - 


12.5 


15.51(8) 15.37(9) 15.39(6) 
15.46(7) 15.36(8) 15.04(13) 


0.417(4) 0.425(5) 
0.435(4) 0.431(8) 


0.81(2) 0.001 
0.84(7) 0.000 



The comparison with the PIMC results for T\ = 1 and 
T 2 = 0.714 is shown in Fig [6] Two system sizes (Vi,^) 
are considered to demonstrate the finite size effects. They 
are found to be negligible and are within the statistical er- 
rors of n q . The agreement with the asymptotic (48 1 can 
be confirmed for the weak (D = 0.1,0.5) and intermedi- 
ate coupling (D — 1.75). For D > 7.5 we observe the 
onset of the slope predicted by (48 1. 
(V > 10 3 ) are required to access smaller ^-values. 



c and the superfluid fraction are given in Tab. VI 



Larger system sizes 
Inde- 
pendently, one can check the reconstructed spectral density 
A(q, uj) in Sec. |VIII A||vT!lC| In the limit q, u> -> 0, for all 
coupling strengths A(q, oj) has a pole at uj = cq, where c co- 
incides with the compressional sound speed. This explains 
a good agreement in Fig. |6^,b. However, for D > 1.75, a 
second excitation branch appears in A(q, lo) (see Fig. [8]). Its 
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spectral weight increases with q and coupling D. This can be 
the reason for the observed systematic deviation in Fig. |6j;,d 
and its onset at smaller q as D is increased. 



B. Collective excitations 



T=l 



T=2 



T=3.3 



As discussed in Sec. IV in the superfluid phase, the single- 



VIII A 



particle (SP) spectrum is coupled to the spectrum of density 
fluctuations. The first term S\(q,uj) in Eq. (10 1 describes 
the scattering of quasiparticles in and out of a condensate, 
with a sharp <5-like peak quasiparticle dispersion expected at 
T <C T c . These sharp features should be present in S(q, lS) 
in the superfluid phase and vanish for T > T c . Hence, it is 
instructive to analyze the T-dependence of S(q, uS) to see any 
difference in both phases. We first discuss some general fea- 
tures observed for different couplings (0.1 < D < 12.5) or 
densities (D ~ s/n), and then go into detail in Sec 
I VIII CI 

The dynamic structure factor reconstructed by the SO 
method is represented in Fig. [7] In the superfluid (T = 1) 
we observe at least two excitation branches with well defined 
dispersions (with the sharp energy resonances in the phonon 
and roton part of the spectrum). With D increasing the onset 
of the second high energy (H) branch systematically shifts to 
smaller q-vectors: for D — 0.1 (D = 12.5) it gets a signifi- 
cant spectral weight at qa > 6 (qa > 3). The if -branch does 
not vanish in the normal phase (T = 2.0,3.3), but is suffi- 
ciently damped or merge with the lower branch. This implies 
that it is closely related to the multiparticle excitations. 

The dynamic structure factor of the upper branch Sh{q> w ) 
(see Fig. [7]) is strongly influenced by the statistical noise in 



the imaginary correlation function G2(q, r) (118b as discussed 
in Appendix |A| The contribution of the high energy features 
is exponentially damped by the factor e _T ". To accurately re- 
solve the form of S(q, ui) at large frequencies requires higher 
accuracy. Moreover, the reconstruction procedure (Sec. |VI|i 



has a tendency to underestimate the half-width of the high- 



frequency structure (Sec. A 2 1. See also a note in Appendix|D] 



In contrast, the low-energy (L) features can be reproduced 
more accurately (Sec. A2i. The L-branch remains in the 



spectrum both at low and high temperatures and shows the 
temperature induced broadening of the half-width of the spec- 
tral peaks characterizing the inverse excitation-lifetime. The 
dispersion remains well defined in the acoustic range of the 
spectrum and gets significant broadening at large wavevec- 
tors (qa > 8). Near the origin (q — 0) the half-width is 
increased again due to the off-resonant excitations by ther- 
mal fluctuations, huj(q) < ksT. For simulated temperatures 
T = 1.0, 2.0 and 3.3 in Fig. |7]this effect is well observed for 
D < 1.75. 

Now we discuss the features which build up due to the cor- 
relation effects. For D > 1.75 and qa > 3 (see Fig.|7^,k,n) 
the L-branch dispersion bends down forming a local maxi- 
mum - a maxon. At larger wavevectors qa E (6, 8) the roton- 
minimum is observed. The critical coupling for the roton for- 
mation, D = 1 .0 — 1.75, is in agreement with the previous 
analysesPSISlAs the present reconstructed spectrum is free 




10 20 30 40 50 60 70 10 20 30 40 50 60 70 10 20 30 40 50 60 70 
Energy, huj/Eo 

Figure 7: Rescaled dynamic structure factor S(q,u>)/ S(q) at dipolar 
coupling D — 0.1,0.5,1.75,7.5,12.5 and three temperatures T. 
The system is superfluid (normal gas) at T — 1 (T = 2, 3.3). The 
discrete wave-numbers, q = 2im/L (n = 0, 1, . . .), are induced by 
the periodic boundary conditions. The dashed line (guide to the eye) 
is the upper bound ui x for the lower (L) excitation branch derived 
from the sum rules (see Sec.lBl. The upper (H) branch lies above. 



of any approximations the roton-depth is found to be lower 
then reported before. In Refill it was found that the upper 
bound estimate uj v (q) (Eq. B20l is at least as good as the cor- 

t (CBF) from RefP in the roton- 



"xC?) (Eq- 

related basis function resu 



region and predicts a deeper roton-minimum for strong cou- 
pling D > 5. The reconstructed dispersion (see Figs. [TUfr 
and [TTfr ) shows that the correct maxon-roton dispersion is 



even lower than 0J x (q). According to the sum rules (|B8|l-( re- 
fcomp) presence of an upper iJ-branch and the increase of its 
spectral weight in S(q, uS) should push the L-branch to lower 
energies and, correspondingly, deepen the roton feature. We 
also do not exclude that the discrepancy with the CBF result 
(T = 0) is a temperature effect. As was shown in Ref.^J for 
superfluid 4 He one observes a softening of the roton mode 
once approaching T c from below, while the peak in S(q, ui) 
shifts to lower frequencies. 
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In 2D dipolar systems the roton-feature is a pure correlation 
effect which cannot be reproduced by the Bogolubov disper- 
sion, (jJ 2 (q) = e 2 + 2no(T)V(q)e q . The Fourier amplitude 
V(q) of the dipole potential 96 is positive at all q-vectors and 
cannot lead to the rotonization of the spectrum a t weak cou- 
pling (low densities), in contrast to 3D geometry ! 32 ! 42 ! 74 ! Even 
at the lowest coupling considered (D = 0.1) the dispersion 
deviates from the Bogolubov result. The basic assumption - a 
small depletion of the condensate, is not satisfied for D > 0.1 
(see Tabjv}. The zero temperature analysis 36 ends up with the 
same conclusion. 

Next, we analyze the splitting of the dispersion curve into 
the L and H -branch. The onset of splitting can be predicted 
based on the f-sum rules (|B7|)-(|B8|l. In the superfluid phase 



both branches are well defined (Fig. |7^,d,g,k,n) and we can 
consider the ansatz 



S(q,ui) = S L {q)8(io - uj L (q)) + S H (q)S(ui 



(51) 



where ojuh) an d Sum define the dispersion and the spectral 
weight for two branches. Substituted in Eqs. (B7 i-(B8 i the 
system of coupled equations can be solved with the result 



ui H (q) = 

S H (q) = 
S L (q) = 



uj F (q) -u L (g) 
l-u L (q)/u) x (qY 
q 2 [l-uj L (q)/uj F (q)} 
2m LO H (q) - u L (q) 
S(q) - S H (q). 



(52) 

(53) 
(54) 



The solutions depend on the dispersion of the L-branch 
ujL(q) assumed to be known from the SO reconstruction. 
Two upper bounds bjp{q) — q 2 /2mS(q,T) and ui x (q) — 
2nS(q,T)/\x{q,T)\ are defined in Appendix |b| and can be 
evaluated via the static structure factor S(q) and the static 
density response function x(o)- The energies of the lower 
branch ojl, due to a slower decay in the imaginary time, can 
be resolved by the SO-reconstruction more accurately than 
ujh and, therefore, are considered as an input. In addition, 
u>l are less damped at large q-vectors compared to the fre- 
quencies of the collective modes. The solution for the up- 
per branch ujh is compared in Figs. |9fr{T0fcfTTfr (solid gray 
line) with the full SO-reconstruction (solid symbols show po- 
sition of the maxima including the half-width). The SO data 
are in good agreement except for the q-vectors when the fre- 
quencies ujl and oj x overlap or the Lf-branch is significantly 
damped. Therefore, such analysis is not applied at T = 2.0 
and 3.3. In summary, the SO-spectrum provides more infor- 
mation (and more complicated structure) than suggested by 
the simple ansatz pi} . New parameters (additional excitation 
branches and their spectral weights) can be added in Eq. ( f5T| , 
however, to resolve them one needs to know the additional 
frequency moments (uj n ). 

The restriction on the spectral weight to be positive, i.e. 
Sii(q) > 0, predicts the first appearance (at a specific q- 
vector) of the high energy branch ujjj in the spectrum. This 
requires cu F (q*) > uj L (q*), see Eq. (53 1. For D = 0.1 
this occurs at q*a sa 5, for D = 1.75(12.5) at q*a w 



1.5(1.0). This comes in agreement with the SO spectrum, see 
Figs. |9frfT0frfTTfc, and confirms the self-consistency of the re- 
constructed S(q,u) with the /-sum rules (B7i-(B9i. The ac- 
curacy in the fulfillment of (B7 i-(B9 1 evaluated from S(q, ui) 
varies in the range 10~ 5 . . . 10 -3 . 



The solutions (52 1-(53 1 can also predict a q-dependence of 



the spectral weights. For D = 0.1, Su(q) is enhanced around 
qa ps 6.36 with u>u being slightly above the recoil energy e q . 
Next, Su(q) is slightly increasing while the intensity of the 
L-branch is decreasing. The SO data confirm this variation, 
see Fig. [7^. Such a behavior is imposed by almost a constant 
value of S(q) in Eq. (54 1 for large momenta. 



The theoretical interpretation of the splitting into two 
branches will be further discussed in Sec. I VIII Al Here we 
mention two possible scenarios, considering as an example 
the dispersion for D = 0.1 (Fig. [7^): I) After hybridization 
at qa s» 6.4 the observed lower branch Wi(g) is the con- 
tinuation of the a single-particle (SP) dispersion or a disper- 
sion of collective density modes. The high energy resonances 
(H -branch) correspond to combinations of two or more quasi- 
particle excitations, II) both dispersions of the SP and collec- 
tive modes for qa > 6.4 continue above (not lower) the free- 
particle dispersion e q . The L-branch appears due to decay pro- 
cesses, then the quasi-particle energy becomes larger than the 
energy of a quasiparticle pair, i.e. u>z(q) > u>z(qi) + u)L,{q2) 
with q = q 1 + q 2 . 



C. Single-particle excitations 

The reconstructed spectral density \A(q, uj)\ is demon- 
strated in Fig. [8] The density at positive A > (q,u>) > and 
negative A < (q,ui) < frequencies satisfies the normaliza- 
tion (|3~8]l-(|39]l and characterizes the excitations that are ei- 
ther internally excited by thermal and quantum fluctuations 
(ui < 0) or excited by the energy transfer (ui > 0). The de- 
cay of the negative amplitude \A<\ with the q-vector is due 



to the normalization (39i and follows the decay of the mo- 



mentum distribution n q = Gi(q,(3) depending both on the 
dipole strength D and temperature T (see Fig. |5J. The weak 
T-dependence of n q observed for D > 1.75 and qa > 1 can 
be directly linked with Ay and A<. Rewritten in the forrrP^I 



dw 
dw 



-1] 1 [A > ( g , W ) + A < (( 7 ,- W )] 



A<(q,u) 



(55) 



the T-dependence enters through the first antisymmetric com- 
ponent. If the spectral densities Ay and A< are antisymmet- 
ric, the second term will dominate being proportional to n q . In 
this case the T-dependence does not enter explicitly and ap- 
pears only indirectly via multi -excitation and damping effects 
contained in A< (q, to). Therefore, the observed symmetry of 
the lowest excitation branches cj>,l, ^<,l for D = 1.75, 12.5 
in Fig. [8] should result in a weak T-dependence of the mo- 
mentum distribution, which is in agreement with n q (T) in 
Fig.|5j;,e. 
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Figure 8: Spectral density \A(q, u)\ at D = 0.1, 0.5, 1.75, 12.5 and 
three temperatures. The system is superfluid at T = 1. The dashed 
line is the upper bound uj x for the ljl -branch in S(q, cj). The spec- 
tral densities are not smoothed and demonstrate the raw-data. The 
frequency resolution compared to Fig.|7]is lowered from w m j n = 0.5 
to Wrma = 1.5 (see Sec.[VT]for details). Hence, the sharpest features 
have the smallest possible width equal lUmin. 



The SP excitation spectrum has a complicated structure and 
many high-frequency harmonics compared to the collective 
density modes. For the strong coupling D = 12.5, to fit the 
imaginary time decay of the Matsubara Green function ( 17 i 
the reconstruction was performed on the enlarged frequency 
interval. Some high-frequency harmonics are out of the scale 
shown in Fig. [8] The detailed interpretation of the SP spec- 
trum and check of the corresponding /-sum rules is out of the 
scope of the present paper and requires further analyses. 

Below we discuss the differences observed in the recon- 
structed spectral densities S(q,u>) and A(q,u) in the super- 
fluid and normal phases in three coupling regimes. Possible 
experimental realizations were listed in Sec. [Ill] 



VIII. DISCUSSION OF THE SPECTRA 

A. Weak coupling 

First, we analyze S(q, u) at D = 0.1. In Fig.|9^-c we plot 
the position of the peaks and their half-width. The black-solid 



symbols denote the main spectral features in S(q,u>). Three 
temperatures are compared. In the superfluid phase (T = 1) 
in the range qa < 8 we observe a sharply peaked single ex- 
citation branch. The dispersion is accurately reproduced by 
the upper bound oj x (q) derived from the /-sum rules, see 
Appendix [B] The half-width of the peaks characterize the 
damping effects. The smallest half-width in Fig. [9] is lim- 
ited by the frequency -resolution used in the SO-reconstruction 
(hAu>/Eo = 1, Sec.fVl]). Except of the region, hw(q) < k B T, 
the half-width stays close to this lower bound for qa < 5. 
For larger q the damping starts to increase systematically (see 
T = 2, 3.3). In the region (qa rj 5) the dispersion curve 
crosses the free-particle branch, e q = q 2 /2m, and broadens. 
This can be interpreted as hybridization and level repulsion 
which occurs whenever two branches cross. In the superfluid 
(T = 1) this broadening is only slightly increased with the 
q-vector and simultaneously we find the splitting into the L- 
and iJ-branches. The L-branch goes well below the recoil 
energy and follows the dispersion of the zero sound (ZS). In- 
terestingly, this feature can be only observed in the superfluid 
phase. Based on the theory of the hybridization of G\{q) and 
x(q) in the presence of a condensate (their poles are repro- 
duced in each function), this branch should be due to the cou- 
pling with the SP spectrum. However, at large q this branch 
is not observed in the spectral function A(q, to). For a further 
test, one needs to evaluate the two-particle Green function and 
check the two-particle spectrum where this branch should get 
a significant spectral weight. 

This behavior is not reproduced in the normal phase, see 
Fig.[9]>c. The damping is systematically increasing with the 
q-vector which is typical for a normal gas/liquid. In addi- 
tion, we observe that the L-branch is now damped and sat- 
urates near a "plateau" well below the ZS -dispersion. This 
shift to lower frequencies is accompanied by the shift of the 
ff-dispersion to higher excitation energies. This result is ex- 



pected from the sum rule (B8 1. We remark that the half-width 
of the -branch can be underestimated by the reconstruction 
(see Appendix|A|). 

What is common for both superfluid and normal phase 
is the linear acoustic phonon dispersion, uj(q,T) = c(T)q, 
which extends to qa « 5 at T = 1.0, qa « 4 at T = 2.0 
and qa s» 3 at T — 3.3. The results for the isothermal 



sound speed c(T) are given in Tab. VI Its non-monotonic In- 
dependence (included in Fig. [9]) can be explained by suppres- 
sion of the density fluctuations in the pre-superfluid regime 
(T < 2) and the corresponding non-monotonic behavior of 
the compressibility kt- The variation of the zero sound 97 
5c = |1 - c(2)/c(l)| evaluated from Tab.|VI 
nounced (5c 



is most pro- 
7%) for weak coupling D — 0.1. At larger 



coupling D the density fluctuations due to formation of a lo- 
cal condensate and superfluid density are strongly suppressed 
by correlation effects (the relative change of the compressibil- 
ity is also reduced). 

To test the hybridization with the SP spectrum in the su- 
perfluid, Fig. |9ji-f shows the spectral function A(q, u). The 
statistical error in evaluation of Gi(q, r) is larger then in 
G2(q, t) (compare Fig. [3] and Fig. EJ, as a result the recon- 
structed SP dispersion curve at low g-vectors (Fig. |9ji) shows 
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Figure 9: (Color online) Dispersion relations for D — 0.1: posi- 
tions of the peaks of S(q, ui) (a-c) and A(q, ui) (d-f) and their half- 
width shown as the error-bars. The full spectral densities S(q, ui) 
and A(q, uj) are shown in Fig.[7jand Fig. [8] Temperature T = 1 is 
below and T = 2.0, 3.3 are above the critical temperature T c (see 
Tab. |VI) , The simulation parameters are specified in Tab. [Il] The 
solid (black) line uj x (q) an d me dashed (blue) line ujf(ci) are the 
upper bounds for the lowest excitation branch derived from the com- 
pressibility and first-frequency sum rules (see Appendix|5J. The dot- 
ted (brown) lines are the free-particle (q 2 /2m) and isothermal sound 
(cq) dispersions (c is taken from Tab. |VI| l. The solid gray line is the 
solution of Eqs. 1 52 i-i 54 1 for the high energy branch ujn(q) plotted 



for the q- vectors where the spectral weight Sn(q) is positive. The 
red symbols denote some additional features in the spectrum which 
are hardly seen in Figs. 17] [8] 



some statistical fluctuations around the ZS dispersion. At 
high temperatures the convergence to the linear dispersion 
(in the range qa < 5) is obscured by the increased damping 
(Fig. |9^,f). However, at larger q the accuracy of the recon- 
struction should be improved as the Matsubara Green func- 
tion is evaluated more accurately with less statistical noise. In 
comparison with the collective excitations, the overall slope 
of the SP dispersion is less influenced by temperature. 

The linear dispersion observed as q — > is in agreement 
with the Hugenholtz-Pines sum rule for the self energies 73 
with the result that the spectrum A(q, uj) (at T = 0) is gap- 
less in the long-wavelength limit. Similar result for Bose liq- 
uids at T 7^ has been worked out by Cheung and GriffinP^ 
Due to the hybridization the density-density response and one- 
particle Green functions in the superfluid share the same sin- 
gularities. At q — > these are the phonon poles as was shown 
by the field- theoretical calculations of Gavoret and Nozieres 46 
with the velocity precisely equal to the ZS speed. Thi s resu lt is 
confirmed for superfluid helium both experimentally^^ anc j 



theoreticallyP^EnEIl This justifies the Landau-Feynman in- 
terpretation of the sound waves as elementary excitations in 
Bose systems. Excited out of a Bose condensate a quasiparti- 
cle excites the collective oscillations in a superfluid. 

In the dipole system, the linear dispersion terminates at 
qa 3 — 5.4 (depending on T) with the splitting and broaden- 
ing. Then the dispersion curve goes slightly above the recoil 
energy e q . Fig. [9^ shows that the same branch (iJ-branch) 
is observed in the dynamic structure factor. This comparison 
confirms, that, indeed, in the superfluid phase both spectra are 
coupled. In the superfluid (Fig.|9^,d) at D = 0.1 both disper- 
sions are linear and coincide with the compressional sound. 
The broad high-energy branch at T > 2 and large wavevec- 
tors (qa > 5.4) is identified with the multiparticle excitations. 
The halfwidth of this dispersion is reduced in the superfluid 
(Fig. [9^) but not drastically. Still the damping is much larger 
compared with the lifetime of the SP-mode in Fig.|9jl. 

Finally, we interpret the L-branch in Fig. [9^ as a continua- 
tion of the ZS mode which terminates (loses its intensity) for 
qa > 9. 



B. Intermediate coupling 



In Sec. IVIIBI we discussed formation of the maxon-roton 
branch for D = 1 ... 1.75 and failure of the Bogolubov 
spectrum. At T = 1 and qa < 3 the dispersion rela- 
tion is accurately reproduced by us x , see Fig. 10 1. Simi- 
lar to D = 0.1 the deviations appear together with the H- 
branch. This branch has a weak g-dependence and spans the 
frequency range ujh ~ 20 ... 30. The resonant frequencies 
can be explained by a combination of low laying excitations, 
U}(qi)+uj(q2) « 10 + 10 R3 20. The g-dependence is in agree- 
ment with the solution ujn(q) from Eq. (52 1. At qa > 8 there 
is an abrupt broadening of the linewidths. For qa > 9 the 
L- and H -branches merge into a single dispersion which goes 
slightly below e q , see Fig.fTOk. The peak position is accurately 



predicted by ui x . A similar shift below e q was experimentally 
observed for liquid 4 He^4 



In the presence of the H -branch (Fig. 10 1) the low-energy 
dispersion has a well defined maxon-roton feature at q E 
(6, 8). I ts depth is lower than found in the previous analy- 
ses j 30 | 36 | -j^jg p resen t approach is not limited by three-particle 
decay processes of CBF theory 3642 and from first principles 
includes the spectral weights and energies of all excitation 
branches, here the L- and H -branch. They are coupled via 



the /-sum rules ( B7 1-( B9 1 as in the simple example ( 52 1-( 54 1. 
Therefore, the under(over)estimation of one of the branches 
has a strong influence on the full spectrum. 

In agreement with Gavoret and NozieresP^in the q, uj — s- 
limit the density response and SP spectra converge to the ZS 



dispersion. The sound speed c is given in Tab. VI At temper- 
ature T = 1 in the SP spectrum, we observe, in addition, a 
second excitation branch with a finite gap ui g « 40 at q = 0, 



see Figs. R5g, 101. The gap value can be explained by the 



excitation of two quasiparticle with the opposite momentum, 
u(q) + oj(-q) m 40 with ui(q) w 20 for qa £ (2, 6). 

We can confirm that the lower dispersions in S(q, uj) and 



17 



A(q, (J) coincide very well up to a maxon region (Fig. 10 i,d) 



T=l 



T=2 



T=3.3 



In the maxon-roton region, both spectra are not much simi- 
lar. Starting from qa ~ 4 the i7-branch in Fig. 10 1 bends up, 



while the SP dispersion in Fig. 101 goes down approaching 
a local minimum near the roton wavevector qa sw 7. Both 
spectra coincide again only for large momenta and follow the 
lu x -dispersion. From these results we can not confirm exis- 
tence of the unified excitation branch in the whole range of 
q-vectors, as was discussed in Sec. |l]and argued in Refill 



Here, we note that the SP spectra in Fig. 101-f are more 
difficult to reconstruct as accurately as S(q, to). The statistical 
noise of the Matsubara Green function G\ (q, r) is by factor 10 
larger than in G\(q, t) (see Fig. [3j, while positions of the en- 
ergy resonances are influenced by the noise level as is shown 
in Appendix |A| For A(q,u) one also needs to accurately re- 
construct the additional resonances at lu < (Fig.lSll. This can 



be a reason why for some q-vectors the spectrum in Fig. 101-f 
does not behave like a continuous dispersion relation. 

Next, we continue to discuss the collective excitations in 
the normal phase, Fig. 10 3,c. The splitting of the dispersion 
curve, previously observed for D = 0.1, here is also repro- 
duced. At T = 2 the roton peak shifts to low frequencies 
(Fig.fToT^). Similar behavior was observed in the normal phase 
of 4 HeJ bl l 62 l where above T c the softening of the roton mode 
was found with the roton energy approaching the zero fre- 
quency. At high temperature (T = 3) the roton-feature gets 
broader and eventually becomes nearly dispersionless and sat- 
urates as a "plateau" (Fig. [10];). Simultaneously the upper 
branch demonstrates a large damping and gradually trans- 
forms into the free-particle dispersion. In conclusion, our 
analyses show that the roton mode, which is very sharp in the 
superfluid phase (T = 1), remains also in the normal phase of 
the dipolar gas. Its lifetime drops significantly when the tem- 
perature is increased from T = 1 to T = 2. Then the damping 
stays practically constant by further increase to T = 3. 

The final form of the dispersion (lower branch) is typi- 
cal for a normal gas/fluid. After the phonon part the dis- 
persion slightly bends down. The self-consistent field pro- 
duced by the density fluctuations cannot support stable short- 
wavelength collective modes. In the normal phase the disper- 
sion relation can be compared with the results of MD simula- 
tions for 2D classical dipoles^J and predictions of the quasi- 
localized charge approximation, 100 after the correct mapping 
of the quantum dipolar coupling D — p 2 m/ e\,h 2 a on the clas- 
sical coupling parameter = p 2 /a^ksT. 



Interestingly, that in the SP spectrum in Fig. 10 1 also re- 
sembles a roton mode but with larger excitation energies, 
compared to the roton in the density fluctuation spectrum. The 
slope of the L-branch is close to ui x (see Fig. [8}. Besides the 
L-branch, there is a well distinguished upper branch. This is 
a common spectral feature for all coupling strengths (Fig. [8}. 
The Lf -branch is very pronounced for D > 1.75 and accumu- 
lates most of the spectral density for u > 0. 




Figure 10: (Color online) Dispersion relations S(q, to) (a-c) and 
A(q, uj) (d-f) for D = 1.75. The half-width of the peaks is shown as 
the error-bars. See Fig.[9]for further details. 



C. Strong coupling and rotonization 

As a third example we consider strong coupling. First, the 
dispersion in the superfluid phase will be discussed. For sim- 
ulated temperatures, at D = 7.5 and 12.5, there is a local 
crystalline ordering with fluctuating orientation. This regime 
is followed by the freezing transition at D = 17(1 ) T232H 
At D = 12.5 the static structure is peaked at qa = 6.65 
and the phonon-maxon-roton feature is well pronounced, see 
Fig.|7]<;,n. Due the 0— sum rule (|B7|i the rotons are the domi- 



nant excitations (the integrated spectral density is proportional 
to S(q)) and the depth of the roton minimum affects the crit- 
ical temperature of the superfluid transition. 30 ^ The full re- 
construction clearly demonstrates the temperature effect, see 
Fig. [7ji. Near the roton-minimum at Er = wl(6.36) ~ 4.9 
the dynamic structure factor S(q,co) shows an asymmetric 
broadening to lower energies. This also applies to D = 7.5 
(Fig. [7jc) with the roton gap Er = wl(6.36) « 6.95. For 
D = 12.5, the simulated temperature T = 1 is close to the 
BKT-temperature of an infinite system (T c — 1.01) and the 
superfluid fraction is reduced to p s = 0.81 (see Tab. VI 1. 
Therefore, we observe some additional damping compared to 
D = 7.5 with p s = 0.95 and T c = 1.22. 

Except the maxon region, the upper bound io x (q) repro- 
duces well the L-branch for qa < 8 including the roton mini- 
mum, see Fig.[TTfr,b,c. The deviations at large q are due to de- 
cay processes. Once the quasiparticle energy exceeds the en- 
ergy of a quasiparticle pair, it becomes unstable and would de- 
cay into this pair. For a strongly correlated system the lowest 
quasiparticle energy corresponds to a roton, Er, and, hence, 
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the dispersion curve should not exceed 2Er. 

Surprisingly, at D = 12.5 the roton minimum is so deep, 
that the energy of the two-roton state 2Er w 10 is much 
lower than the maxon energy. As a result we observe this 
state symmetrically on both sides of the roton minimum at 
qa s» 6.36 (see horizontal dashed lines in Fig.fPTfr). A similar 
effect has been discussed for 4 He at high pressures (above 18 
bars)P^For D = 7.5 (see Fig . |7jc) the double roton-feature 
at 2Er ps 14 is also reproduced. The two-roton state, first 
predicted by Pitaevskiipil is well known experimentally for 
superfiuid 4 He and was partially discussed in Sec.|l] Its exper- 
imental verification was complicated, since the sharp energy 
resonance merges with a broad multi-excitation background 
and the measurements required a high instrumental resolu- 
tion!^ In contrast, for 2D dipolar gases with D > 7.5 there 
is no such complication, since both components are well sep- 
arated (see Figs.fTTfr and[7]<:,n). 
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Next, we discuss the maxon region, qa € (2, 4) (Fig. 11 i). 
Here, we observe a splitting of the dispersion. The lowest 
branch is the two-roton state. Its energy is well below U) x (q) 
(solid black curve) which corresponds to a prediction of a con- 
tinuous dispersion which neglects decay processes. It pro- 
vides an accurate prediction only up to the maxon. Here, the 
upper branch starts which continues with few oscillations (see 
Fig. [7ji). In the roton region, q £ (5,7), its energy is in the 
range u s» 20—25 and there is a large gap to the low frequency 
roton state. At large momenta the ii-branch merges with the 
strongly damped dispersion of multiparticle excitations near 
the recoil energy. Below we give a possible interpretation of 
the observed behavior. 

In the interpretation proposed by Glyde and Griffin 4 ^ the 
hybridization of the dispersion near a maxon is due to the 
crossing of two different branches, i.e the acoustic phonons 
which dominate at low q and the maxon-roton mode. This 
cross-over behavior should be characterized by a double peak 
structure. Indeed, in the simulations we observe a double peak 
structure, but its origin is different. It comes from the decay 
of quasiparticles into pairs of lower energy (here the two-roton 
state). If this decay is not allowed by kinematics then we al- 
ways observe a continuous dispersion relation (see Fig. [9^-c 
and 10 i-c), in agreement with Ref.^The dispersion does not 




lose its spectral intensity or increases its halfwidth as it would 
be for the case for hybridization of two distinct branches. Fur- 
ther analyses for the case when 2Er > E maxon will be helpfull 
to finally clarify wether the maxon for strongly coupled sys- 
tems is a transition region between phonons and rotons or a 
part of a continuous dispersion relation. 

Another important issue concerns the origin of the roton 
minimum. For D = 1.75 we have observed a weak rotoniza- 
tion of the spectrum both in the superfiuid and normal phases. 
Its existence is independent on the presence of a condensate. 
For strong coupling, D = 12.5, the maxon-roton feature is 
more pronounced including the normal phase (Fig. [TTp,c). 
Temperature, practically, does not influence the slope of the 
low-frequency roton-two roton state, but has a noticeable ef- 
fect on the damping characterized by the half-width. We can 
conclude, that the roton minimum is not related with the cou- 
pling to the SP excitation spectrum, being an intrinsic col- 
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Figure 11: (Color online) Dispersion relations S(q,to) (a-c) and 
A(q,u)) (d-f) for D — 12.5. The half-width of the peaks is shown 
by error-bars. See Fig. [9]for further details. In d)-f) the black color 
denotes the uppermost branch which is omitted from the discussion 
in the text. 



lective density response mode. Its formation is solely due 
to th e sho rt-range correlations which supports earlier sugges- 
tion d 70 ! 71 ! and is in agreement with the analysis of 2D dipolar 
system in the classical limitP 

The next question to be addressed is whether the roton 
feature can be also found in the SP spectrum A(q,u). The 
peak positions of A(q, to) and their half-widths are shown in 
Fig.[TTjl-f. We do not discuss the uppermost branch. In con- 
trast to S(q,Lu), the high frequency branches carry the main 
part of the spectral density both in the superfiuid and normal 
phases. They also dominate over the i-branch in the roton- 
region, see Figs. [SJi-p. Two lower branches are well sepa- 
rated at T = 1, see Fig. [PTft . Temperature systematically 
reduces the energy of the upper branch (H) and flattens the 
lower branch (L). 

In the superfiuid phase (Fig. [TTJl) the L-branch has a 
roton-like dispersion, resembling the density response roton 
in S(q, uj). Similar to the density -roton, the left roton-branch 
does not reach the maxon maximum (shown by w x ), but for 
qa e (2, 4) decays into quasiparticles with lower energy. 

The T-dependence of the roton in A(q, w) is different from 
the density-roton in S(q,u>). Both in the superfiuid and nor- 
mal phases it stays very near or above u x . There is no con- 
tradiction as uj x is the upper bound for the collective den- 
sity modes. The dispersion becomes almost flat at T = 3.3 
(Fig. [TTf ) and the roton mode completely vanishes. In con- 
trast, the density-roton mode and two-roton states do not 
change much except for the broadening. In the long wave- 
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length the dispersion converges to the ZS. The lowest branch 
is gapless (w, q —> 0) in agreement with the Hugenholtz-Pines 
sum rule. 

On the right side (qa > 7) the roton-branch follows the 
dispersion predicted by ui x for the density modes. Simulta- 
neously, we observe that there are no excitations at the free- 
particle dispersion e q . This is in contrast, with the behavior 
in the normal phase (Fig. [TTJ,f). Now, the density of the SP 
excitations is centered around the free -particle branch e q . The 
dispersion is broadened significantly compared with the one 
in the superfiuid phase. From this result we can argue that 
the hybridization scenario, discussed in Sec. [I] is realized in 
our 2D system of bosonic dipoles. The poles of the dynamic 
response function x(l> w ) predicted based on the sum rules 
result co x (q) (see Sec. [Bji can be directly observed in the el- 
ementary excitation spectrum A(q, ui) but only in the super- 
fluid phase, when this is allowed by the hybridization. On the 
other side, we do not see exactly the same poles in S(q, w), as 
the collective spectrum is strongly renormalized (perturbed) 
by the quasiparticle decay processes. They shift the positions 
of the energy resonances and lead to a slightly different dis- 
persion relation. This can be a possible explanation for our 
numerical results. 

Finally, we note that many additional modes are observed in 
A(q, lS) which, probably, are intrinsic to the single-particle ex- 
citations. Their interpretation requires further analysis which 
goes beyond the present work. 



IX. SUMMARY 

Using the path integral Monte Carlo technique in the grand 
canonical ensemble we evaluated the imaginary time density 
response and the one-particle Matsubara Green's function. 
The stochastic optimization method was used to accurately 
reconstruct the dynamic structure factor S(q,uj) and the SP 
spectral density A(q, uS) in a wide range of wavevectors. 

As a result, a first-principle treatment of the dynamic struc- 
ture factor was presented both for weakly and strongly in- 
teracting Bose-condensed systems. The temperature-density 
dependence of the phonon-maxon-roton dispersion, its damp- 
ing and hybridization was resolved by the simulations. The 
dispersion relations were compared vs. the upper bound de- 
rived from the frequency sum rules. The introduced approach 
allows to discuss the physical origin of the roton minimum 
and is a useful benchmark for theoretical studies and interpre- 
tations of experiments on quantum liquids/gases in different 
physical realizations. As an example, a 2D dipolar Bose gas 
at weak, intermediate and strong coupling has been analyzed 
in detail. 

Our goal was to demonstrate that the single-particle and 
density-response spectra are coupled and share common poles 
due to h ybridiz ation in a superfiuid, as was predicted the- 
oretically! 46 ! 49 ! 50 ! 73 ! Once dipolar bosons become superfiuid 
(T < T c ), we observe a similar dispersion relation in both 
spectral densities A(q,u>) and S(q,u>). However, we cannot 
confirm that this dispersion can be observed at all wavevectors 
as a unified excitation branch 68 when the coupling is strong. 



In this case, the collective spectrum S(q, us) is strongly renor- 
malized by the quasiparticle decay processes. This, according 
to the frequency sum rules, should necessarily lead to a shift 
of the S(q, cj)-maxima from the position of the original en- 
ergy resonances. One prominent example is the decay into 
the two-roton state observed in the strongly correlated Bose 
gas. 

Further, we found that the two-roton state exists indepen- 
dently of the presence of a condensate and, therefore, con- 
clude that the roton mode has a classical origin due to short- 
range particle correlations. Its observation in classical li quids 
is, probably, prevented by the softening and overdamping! 6 -^ 

For quantum dipolar liquids at high coupling/density (here 
D > 7.5) the reconstructed S(q,ui) predicts a pronounced 
roton minimum. The sharp resonant roton feature is trans- 
fered from a superfiuid to a normal phase but with a signifi- 
cant damping. Our reconstruction of the SP spectrum A(q, lS) 
in a superfiuid also predicts the roton-like feature which, how- 
ever, vanishes in the normal phase. Here it is substituted by 
the excitation mode peaked around the recoil energy. 

The common features of the density fluctuation spectrum 
S(q, uj) are the following. For large momenta (qa > 8) we 
observe a strongly damped mode near the recoil energy both 
in the superfiuid and normal phases. As expected, the phonon 
dispersion at small wavevectors exists independently of the 
presence of spatial coherence in the one-particle density ma- 
trix. 

These interesting aspects of the excitation spectrum near 
T c may be addressed in the near future in the ongoing exper- 
iments on dipolar systems. 3 4 The analyzed dipolar coupling 
strengths can be realized in 2D systems of magnetic atoms, 
polar molecules and indirect excitons as discussed in Sec. [HI] 
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Appendix A: Model spectral densities: accuracy of the 
reconstruction 

To demonstrate accuracy of the reconstruction procedure, 
below we consider a test for several spectral densities S(ui). 
The goal is to analyze how a level of statistical noise, present 
in the imaginary correlation function, influences positions of 
energy resonances and their half-width. To characterize devi- 
ations, in addition, we evaluate the frequency moments (ui n ) 
to check fulfillment of the sum rules ( |B7| l-( |B9| >. In our model, 
they are substituted by the exact integral properties of S(u>). 
Here, we analyze a dynamic structure model and assume that 
the detailed balance holds, S(—u) = e^^S^). 

Further, for convenience, S(u>) is substituted by a similar 



20 



distribution S°(lu) discretized into rectangles of width Aw 

^MLe^-iA^+iA^ = S(ui), Ui = (i- ^)Aw. 

(Al) 

By discretization evaluation of the correlation function can be 
easily performed as 



A' 



El „,.+ ., Aut 
-e~ Uit smh 
t 2 

t=T,/3-T 



(A2) 



As a second benefit, the integral properties (frequency mo- 
ments) in the basis of rectangles {hi = S(uJi), Wi = Aw, a = 
uji} can be written down explicitly and evaluated via 



A 



i=l 
N 



Wi + ^-e P Ci sinh 

(3 2 



(A3) 



2hi 



: hjWjCj -pe 0Ci 

»=i " 

~2~' 



1 



sinh 



Wij3 



cosh 



(A4) 



iV 



<o=E^ ln — 



Cj + Wj/2 
Ci - Wi/2 



- hi [Ei(-p(a + wi/2)) - m(-/3(a - te</2))] , 

(A5) 

with Ei(x) = f* dte t /t. To avoid the divergence in 
Eq. (A5 i when (c\ — W\/2) — 0, we set 5° to be non-zero 
for u> £ [6, ujco + S] (with 5 = 10~ 3 and the cut-off frequency 
cjco = 100) and choose, correspondingly, u>i = (i— |)Au+<5. 

Eqs. ( |A3[ )-( |A"5) l can be used to evaluate the frequency mo- 
ments (/-moments) of the ensemble average (25 1 by applying 
to each term in the ensemble. This result was used to check the 
sum rules of the reconstructed S(q, oj) for 2D dipolar bosons. 
However, as we discuss below, an accurate fulfillment some 
of the sum rules [relative error 10 -5 . . . 10 -3 ] does not neces- 
sarily guarantee a correct reconstruction of the shape of S(ui). 

In our tests S°(oj) is specified by a linear combination of 
M Gaussians 



S(f)\u 



1 M 

3=1 



Pj 



2na 



3 = 1 



(A6) 



which vary in the peak position Qj, the variance <jj and the 
weight factors pj . Similar functions has been used in Ref! 3 ^ 
to test a GIFT reconstruction procedure (the genetic inversion 
via falsification of theories). We refer to RefPSfor the comple- 
mentary tests and discussion of limitations by reconstruction 
from imaginary times. 

To simulate the statistical noise inherent to QMC simula- 
tions, the imaginary time correlation function is perturbed by 
a random noise v 



G( n ) = G°{t z ) + 7 ■ v5G{ n ), n = i8r 



(A7) 



sampled from the normal Gaussian distribution and multiplied 
by the statistical error 5G(ri) at imaginary time t,-. The time 
dependence JG(Tj) was taken from a typical PIMC run (e.g. 
a dipolar gas at I? = 1.75 and (5 = 1) and rescaled by factor 
7 (1/5 < 7 < 5) to model the situation when the recon- 
struction is performed for a larger or smaller MC sample size. 
Following typical PIMC data presented in Fig. [2] as setup pa- 
rameters we used: the inverse temperature j3 = 1, the time 
step St = 1/100 (n = iSr G [0, 0/2] and i = ... 50) and 
the statistical error SG(n) of the order 1 0~ 5 . . . 10~ 3 [with 
G(t = 0) ~ 1 due to normalization in Eq. ( A6 1]. 



How particular spectral features can be downgraded by sta- 
tistical errors is clarified by few examples. 



1. Spectral density with double peak 

This example is common for experimental and theoreti- 
cal studies of superfluid helium when two distinct excitation 
branches can be resolved in the low temperature (T < T c ) 
spectrum. The reference spectral density S°(ui) is specified 
by the following parameters: discretization width in Eq. ( |A1[ ) 
Auj = 0.5, two Gaussians with f2j = 10, f2a = 25, <t\ = 0.5, 
02 = 2.0, pi = 1 and p2 = 3. Parameters used in the re- 
construction [see Eq. ( |2"9| ) in Sec. |VI| : the cut-off frequency 
ljco = 100, the allowed width of rectangles w; > 0.5 (i = 
1, N) in the basis of size N = 80, and M = 200 - 400 solu- 



tions in the ensemble average (25 1. 

Fig. [T2fr-f demonstrates the results. In the error-free case 
(SG(Ti) = 0, first column) the reconstructed spectrum (solid 
red line) quite accurately reproduces the position and half- 
width of the first peak in S°(uj) (step-like curve) even on the 



logarithmic scale, see Fig. 12 1. The center of the second peak 



is accurately captured but the shape is perturbed. The recon- 
struction overestimates the spectral density at the peak center 
and the wings. The corresponding relative deviation from the 
reference correlation function G°(t) is shown in Fig. 12 c and 
is of the order of 10~ 6 . . . 10~ 5 . Here we show the deviations 
of the correlation functions G n (Ti) corresponding to differ- 
ence spectral densities S n (u)) used in the ensemble average 
[only 25 samples are shown from typical n = 200 . . . 400]. 
This is the best result we can achieve with the stochastic op- 
timization method (Sec. |VI| i using as an input the error- free 
correlation function G°(r). 

An increase of the basis size in Eq. ( A2 1 to N = 160 does 
not bring any substantial improvement. The reconstruction re- 
sult is mainly determined by the level of statistical noise. We 



also tried to fix the /-moments ( A3 1-( A5 1 to their reference 



values 



/ s° 



(n = —1,0, 1). They were included in (20 1 via 



a set of Lagrange multipliers and subjected to minimization. 
However, the spectrum shape was downgraded by appearance 
of systematically shifted peaks which lead to an exact fulfill- 
ment of the sum rules at a price of a larger deviation from the 
reference function G°(r) along the imaginary time. Surpris- 
ingly, even if the moments are not fixed, their relative error 
S (uj n ) = |1 — ( w ™),s( w ) / ( wn )so( w ) I i s small, see lower pan- 
els of Fig. [12] 



By turning on the noise in the input data ( A7 1 the peak po 
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sitions are slightly perturbed, see Fig. 12 3,c,e,f. The accuracy 
in the peak position is acceptable even in the worst case when 
the error <5G(r i )/G(Ti) is in the range of 2% for r — > 0.5, 
Fig. [I2},m. However, the half-width of the second peak is 
underestimated and should not be taken seriously in a recon- 
struction from QMC data with similar error level. 



The advantage to use the ensemble average ( 25 1 instead of 
individual solutions S n (uj) becomes evident from the panels 
k-m. While one solution G„(t) has a large error at a partic- 
ular imaginary time, it is reduced by the averaging with other 
solutions. 



2. Bimodal distribution with overlap 

In a second example we consider the bimodal distribu- 
tion consisting of two overlapping Gaussians. A sharp peak 
models a quasi-particle resonance on a broad multi-excitation 
background. The reference spectral densities in Fig. 13 and 
Fig.[l4]are defined, correspondingly, by [fli = 10, = 25, 
(Ji = 0.5, a 2 = 8.0, pi = 1, andp 2 = 3] and [Hi — 10, £1 2 — 
12, (Ti = 0.5, <72 = 6.0, pi — 1, and p2 = 1]. 



As in the first example, the reconstruction in the error-free 
case (first column) accurately reproduces the shape and the 



decay at large frequencies on the log-scale (see Fig. 131 and 
Fig. \l4\ l). The deviation along imaginary time (Fig. |T3jc, 
Fig. 14 1) are within the ran ge 1 0~ 8 . . . 10~ 5 . The errors in the 
/-moments of S(oj) in Fig.jBJi are 5 (w n ) = 1.3 • 10~ 6 , 2.4 • 

,0,1. 



S(lu) and the discretization parameter Auj in Eq. (Al 
step-like shape of the reference distribution S° can not be re- 
produced in the smoothed ensemble average result. This re- 
sults in different /-moments even if the shape of both densities 
is close. 

Second and third columns in Figs. |13|14| show the recon- 
struction with the statistical errors in G(r) (cf. the error bars 
SG(ri) in Fig. 13 i,i and the solid black line 5G(Ti)/G(Ti) in 
panels l,m). 

In Fig. [T3p,c we see a trend similar to the double peak 
structure, cf. Fig. [T2"p,c. There is no appreciable difference 
in the position and the half-width of the first peak, but the 
noise shifts the second peak and results in a systematic un- 
derestimation of its half-width. We can conclude, that with 
similar relative errors 10~ 3 . . . 10~ 2 in the imaginary correla- 
tion function from QMC, a spectral shape of a high-frequency 
branch should not be taken too seriously as it is obscured by 
statistical noise. Neverthe less, the relative error to the refer- 
ence /-moments [see Tab. |VII[ is within the range 5 (u n ) — 
10~ 5 . . . 5-10~ 4 (n = —1, 0, 1) and is surprisingly small. This 
demonstrates a difficulty to qualify a reconstructed spectrum 
based only on fulfillment of few sum rules. We return to this 
point again in the summary. 

Fig.[l4]shows an example when a sharp resonance is placed 
on a broad distribution which extends to low frequencies. 
Again in the error-free case the reconstructed density is re- 
liable including the log-scale. 

It becomes difficult to reconstruct this sample when a 
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Figure 12: (Color online) (a-f) Reconstructed spectral density S(lj) 
(solid red line) compared with the reference S°(ui) (step-like black 
curve) consisting of two energy resonances [(d-f) on the logarithmic 
scale]. Three columns group the results for three cases of the rela- 
tive statistical error 8G(ji)/G (t%) shown by a solid (black) line in 
(1-m). First column (a,d,g,k) is the reconstruction from the error-free 
[8G(t) = 0] correlation function G°(t) shown by dashed (black) 
line in (g-i). The bold symbols with error bars specify the data 
[G(n) ± 5G(n), i — . . . 50] used in the error-biased reconstruc- 
tion (second and third columns). Different solutions G„(t) obtained 
by the stochastic optimization are shown by solid (red) lines but are 
practically indistinguishable on the present scale. Each solution cor- 
responds to a spectral density S n used in the ensemble average (??) 
Only part, r g [0.40,0.50], of the full imaginary time interval 
[0, (3 = 1] is shown, (k-m) Solid (gray) lines show the relative devia- 
tions 1 1 — G n (Ti)/G° (tj ) I from the exact correlation function. Only 
25 from n ~ 200 . . . 400 samples are shown. Note, that the devia- 
tions (panels 1-m) are within the statistical error SG(Ti)/G(ri) [solid 
(black) line]. Horizontal line 10~~ 3 is guide for the eye. Three lower 
panels show the relative error of the ensemble average frequency mo- 
ments (|A3|i-((A5]i, i.e. 5 (u n 



*> = |1 



K>sm 



where 



S(lj) is shown in panels (a-c) and denoted on the x-axis as (a), (b) 
and (c). See text for further details and used parameters. 



tiny noise is introduced. The relative error (cf. Fig. 14 ) 
is by two orders of magnitude smaller (3 • 10~ 4 ) then in 
Fig.[T3}ri, though it sufficient to obscure the sharp resonance, 
cf. Fig. [Tffi , and two distributions merge into one. However, 
the reconstruction does reproduce the half-width of the back- 
ground density and its asymptotic decay with frequency (cf. 
panel f). The errors in the /-moments for the reconstruction 
in Fig. 14 i,b,c are within the range 5 (w" ) = 10 -5 . . . 4 • 10~ 3 



22 



0.25 r 

0.2 ■ 
0.15 

0.1 
0.05 ■ 
■ 



-1 



" d) A 










10 20 30 40 50 60 70 10 20 30 40 50 60 70 10 20 30 40 50 60 70 



20 40 60 80 100 20 40 60 80 100 20 40 60 80 100 




0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 



Figure 13: (Color online) (a-f) Reconstructed spectral density S(lj) 
(solid red line) compared with S°(uj) (step-like black curve), (d-f) 
The same on the logarithmic scale. For further details see text and 
caption of Fig. 12 For the relative errors 5 (a/ 1 ) see text. 
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Figure 14: (Color online) (a-f) Reconstructed spectral density S(lj) 
(solid red line) compared with S° (cj) (step-like black curve), (d-f) 
The same on the logarithmic scale. For further details see text and 
caption of Fig. [12] The plots of correlation functions G n (Ti) are 
omitted. For the relative errors S (ui n ) see text. 



(n 



-1,0,1). 



3. Multipeak structure 

In the last example we consider three Gaussians. The ref- 
erence distribution is specified by the parameters [fii = 18, 
£l 2 = 30, and 3 = 60], [a x = 0.5, a 2 = 6, and cr 3 = 10], 
and [pi = 1, p2 — 3, and p% = 3]. The low frequency peak 
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Figure 15: (Color online) (a-f) Reconstructed spectral density S(u) 
(solid red line) compared with S°(uj) (steplike black curve), (d-f) 
The same on the logarithmic scale. For further details see text and 
caption of Fig.|12| 



is sharp and contains 1/7 [cf. Pi/^Pi] of the full spectral 
weight. 



The reconstruction in the error-free case (cf. Fig. 15 i,d,g,k) 
is surprisingly good. All spectral features, cf. the resonance 
and the broad "continuum" are well reproduced. The relative 
errors in the /-moments are S (u n ) ~ 10 -6 , 10 _9 ,3 • 10~ 5 , 
correspondingly, for n = —1,0,1. The deviations along 
imaginary time are within 10~ 8 . . . 10 -5 (cf. Fig. 15 1). 



In second and third columns of Fig. 15 we study effect of 
the noise. In both cases, the relative error for short (r < 0.1) 
and long (r — > 0.5) imaginary times differs by orders of mag- 
nitude and reach 7% and 30% at r = 0.5 (cf. Fig.[B},m). Still 
the properties of the reference spectral density are not com- 
pletely lost. It is possible qualitatively to reproduce the num- 



ber of peaks and their half-width in Fig. 15 3,e and the remi- 
niscence of the first peak, the overall extension over a broad 
frequency interval and the asymptotic decay in Fig. 15:,f. 



Similar to the previous tests, the high-frequency peaks can be 
shifted from their correct position. The shift does not exceed 
17%. 
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Table VII: Frequency moments of the spectral density S°(uj) 
[Eq. ( |A6) discussed in Sec. A 1 | A 2 | A"3] Models 1-4 corresponds 



to S (w) shown in Figs. 12 ■ 151 Used parameters: the inverse tem- 



perature /? = 1.0, the frequency range ui G [5, 100 + S] (8 = 10 
and the discretization Au = 0.5 [Eq. | |A1} ]. 



5 ) 



Model 



(a; 1 ) 



(a;" 1 ) 



1.33334863 28.5554047 0.0729044683 
1.33299124 28.5563603 0.0790881553 
1.32197070 15.5863431 0.142451835 
1.16665053 48.0955388 0.0351857739 



Several first frequency moments can uniquely specify only a 
short time asymptotic of G, and, in general, all moments are 
required to reconstruct G at the full time scale. Vise versa, an 
exact fit to G(t) at all r necessarily means that all frequency 
moments are uniquely defined. However, if there is some un- 
certainty, i.e. G(t)+5G(t), it can be fulfilled by different ex- 
pansions in (oj n ), where the successive terms mutually com- 



pensate, due to the sign-alternating prefactor in (AlOi. Now 
the shape of S(uj) is not unique. 



Appendix B: Sum rules upper bound for the lower excitation 
branch 



4. Summary 

We conclude, from the performed tests, that the reconstruc- 
tion reliability crucially depends on the relative noise in the in- 
put correlation function, 5G(Ti)/G(Ti). This quantity can be 
used to qualify reconstructed spectra, in particular, presented 



in Sec. VII B based on the downgrade of the spectral features 
analyzed in this section. In general, the relative noise in the 
density-density correlation function increases with r and the 
q-vector due to a faster decay of G(q, r) (see Fig. 4a). As a 
result the reconstructed S(q, uj) in Fig. 7 at large wavevectors 
becomes less accurate compared to the low momentum part 
of the spectrum. 

When only two well separated energy resonances are re- 
constructed, the peak positions remain stable against the noise 



(cf. Fig. 12 1. A noticeable effect is observed in the half-widths 



of the peaks. 

The reconstruction results are more sensitive to the noise 



level for the overlapping bimodal (Figs. 13 and 14 1 and the 
multipeak distribution (Fig. [15} . As expected, details on how 
one individual feature (a peak position and a half-width) is re- 
constructed have a strong influence on the rest of the spectral 
density. Even if the low frequency moments are fulfilled with 
high accuracy 10~ 6 — 1CP 4 , this is not sufficient to guaran- 
tee a reliable reconstruction. A large class of solutions G n (t) 
can fit within the error-bars the correlation function ( |A7) . The 
broader is the reference distribution S° in the frequency space, 
more "degrees of freedom" are allowed to be used in the re- 
construction. They can be distributed in many different ways 
to satisfy the sum rules. As a result the ensemble average can 
smooth specific features of the reference spectrum. 

Next, consider a Taylor expansion of the correlation func- 
tion 



G(r) 



^ r n dG^ 

n=0 



dwe-™5(w), (A8) 



from which follows the expansion in the frequency moments 

(A9) 



<9G (n) 

<« B > = (-i) B -^s-i~o. 



G{r) = £ 



(A10) 



The upper bound for the lower excitation branch with a neg- 
ligible damping can be derived using the decomposition of the 
density-density response function (imaginary part) 



%[x(q, w )] = xi(q, w ) + X2(q, w). 



(Bl) 



It is assumed that the <5-peak quasiparticle contribution with 
the energy u q > 

Xi{q,io) = -irA{q)[6{uj - u q ) + 5(u> + uj q )] (B2) 

is not overlapping with the multi-particle excitations at higher 
energies 



X2(g,w) 



7^0, uj > w q , 
0, < uj < uj„ 



(B3) 



We note that Xi(2) is °dd in ui. The fluctuation dissipation 
theorem and the detailed balance are written as 

9f[x(ff,w)] = -7r[l-e-^]5(g,w), (B4) 
S(q,-uj)=e-P"S(q,u>). (B5) 

The dynamic structure satisfies the frequency moments sum 
rules 



dui ui n S(q, uj) 



(B6) 



In the f ollowing, the 0-sum rule, the longitudinal /—sum 
rul er 9 ! 40 1 and the compressibility sum rule (all valid at any T) 
will be employed 



=S(q,T), 
(uj) = q 2 /2m, 

(u- l ) = ~\x{q,T)\. 



(B7) 
(B8) 

(B9) 



These constraints can be used to characterize the quality of 
experimental/theoretical spectra, i.e. check how well they sat- 
isfy the sum rules. 

The finite-temperature generalization of the zero tempera- 
ture Bijl-Feynman upper-bound for lower excitation branch^! 



n=0 



w g < uj° F (q), uj° F (q) = 



2mS{q,T = 0) 



(B10) 
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can be directly derived from the /-sum rule ( |B8| l. The contri- 
bution of the negative frequencies, S(q, —uj), is important at 
high temperatures, typically, in the normal phase. By substi- 



tution of ( B 1 1-( B5 i in ( B7 1-( B8 i we end up with 



S(q) = A(q)cath^ + S Xa , 



U X2 

2m 



duj (3uj 
X2(9,w)coth— . 

— 7T Z 



A(q) uj q 



du) 



X2(q,u)u. 



According to ( |B3| l, we can write the inequality 



— X2 (q, w) <*> > u q tanh ■ S X2 . 



(Bll) 

(B12) 
(B13) 

(B14) 



The combination of ( |B1 l| i,( |B"T3"j ),( |B14| i leads to the upper 
bound 



Wr?t anh^5(q,T) < 



2 v» > - 2m 
where the Feynman frequency ujp is defined by 

(Mq)\ 



or oj q < WF(q), (B15) 



uj F {q) 



2mS(q,T) COth \ 2 



(B16) 



An improved upper bound can be derived from the 
compressibility-sum rule. It was discussed in Ref. 41 and ap- 
plied to 2D dipolar gases in Ref. 30 By substitution of ( BT|b~5 i 
in dB9b we obtain 



\ X (q,T)\ _ + f°° doj X2 {q^) 



2n 



du x 2 (g,w) 



-7T CO 

< w -i tanh ^ • S y 



-7T W 



(B17) 
(B18) 



where we have used that u> 1 tanh[/8a;/2] decays with uj. This 
leads to the upper bound 



with 



0J x (q) — 2nS(q 1 T) tanh 



Puj x {q) 



/\X(Q,T)\. 



(B19) 



(B20) 



The involved static response function can be directly es- 
timated from the imaginary time density correlation func- 
tion (TT~8l 



x(q, T )/ n = - G 2 (q,T)dr. 
Jo 



(B21) 



For a dispersion consisting of a single sharp excitation branch 
the upper bound uj x becomes almost exact, e.g. D = 0.1 and 
Fig. For other couplings, the phonon-maxon and roton 
features are also accurately reproduced (Fig. [9^, [TUfc , [TTfr). 
The reason is the damping of high-energy contributions by 



the factor 1/to. This significantly improves the Feynman up- 
per bound ( |B16| l, where the spectral weight of high frequen- 
cies scales as u>. This leads to the general relation, 0J q (T) < 
0J x (q, T) < uiF(q, T), valid in a superfluid phase for the lower 
excitation branch. 

In the long-wavelength limit, we expect that both up- 
per bounds converge to the isothermal sound dispersion, i.e. 



,(q) « 0JF(q) ~ cq. Then the product of ( B16[ ) and (B20l 



results in the compressibility sum rule 

KmJx(<?,T)| 



mc 2 (T) 



(B22) 



Similar, for fiui < 1 (or cq < k B T), from Eq. ( |B20| i we 
^et the limiting value of the static structure factor 



lim S(q,T) 

q^>0 



lim 

<j->0 



\X(Q,T)\ 

n(3 



k B T 
mc 2 (T) 



(B23) 



The offset from zero increases with temperature and decreases 
with the coupling strength/density (the sound speed increases 
with A seeTab.lVn. 



Appendix C: Bose condensate in 2D excitonic system 

Consider indirect excitons with the dipole length dx = 
20 nm at density n in a circular trap with the diameter 
d = 20 /.im realized in two coupled GaAs QWs structure 
160/40/160 A, as in typical realizations. 8 10 The dipole cou- 
pling can be expressed in the terms of relevant semiconductor 
parameters as 



D 



d 2 



mx 
m* 



a B 



h 2 e 



(CI) 



where for GaAs we take e = 12.58, the in-plane elec- 
tron(hole) mass Ti*/^ = 0.0667(0. 112)m e , the Bohr radius 
a B = 10 nm, the exciton mass m^/m* = 1 + m* h /m* e — 
2.68. The density n 1(2) = 0.22(2.69) • 10 10 cm- 2 will corre- 
spond to the inter-exciton separation au2) = 214(61) nm and 
the coupling Dx(2) = 0.5(1.75). The critical temperature can 
be taken from RefP, T c - 1AE with E = ^[mxa 2 ]- 1 = 
0.0086(0.106) K for the considered densities. Further, we 
take into account that in the unpolarized system the spin- 
degeneracy factor (g = 4) lowers the effective density of 
particles, h = n/g, of the same spin projection which can 
undergo condensation, and, correspondingly, due to the BKT- 
scaling f c - n, we estimate f c = T c /g = 0.003(0.037) K. 
The system size equals the trap diameter L = d = 20/im. 
This corresponds to ^(2) ~ nL 2 ~ 9- 10 3 (1.1 • 10 5 ) excitons 
in the trap. For T\ < T c , the lower bound for the condensate 
wiU be given by n (Ti) > G\ /2 {T{) = 0.38(0.16). For this 
estimation we used the critical exponents from Tab. [II] and the 
reduced system size L/au2) — 94(328). 

To observe quantum degeneracy effects, temperature 
should not be too low (T < T c ). The momentum distribu- 
tion in the pre-condensation regime also has a peculiar shape 
at small momentum (see Fig.[5ji,b,c) for T = 2E /g (T ~ T c ) 
if different spin states are included. 
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Appendix D: High-energy branch 

A simple alternative to reconstruct the dynamic structure 
factor at high energies, S(q, oj)\ui>u*< could be a fit to a 
suggested analytical form, e.g. to a Lorentzian or Gaussian 
g — /a q ^ w j tn tne gQ p rocec i ure operating with the opti- 
mization parameters (uj q , a q ). The reconstructed spectrum in 
Fig.[7]at high T allows such treatment. Then, the SO method 
would reduce to the regularization method, like the ME. How- 
ever, there is a priori assumption that the high-energy branch 



is smooth and broad (in frequency), and there is no overlap 
with (5-like resonances. The single-particle spectral density 
A(q, ui) in Fig. [^demonstrates the potential difficulty with this 
approach when we suggest too simple analytical form con- 
trolled by only few optimization parameters. Therefore, the 
spectral densities are reconstructed by the linear combination 
of rectangles, see Eq. ( |30| l, in the full frequency range with- 
out a priori assumption. This leads to the observed statistical 
noise in the reconstructed densities at high frequencies. 
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